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Introduction 

Approximately  2.5  million  Americans  live  with  epilepsy  and  epilepsy-related  deficits  today,  more 
than  disabled  by  Parkinson  disease  or  brain  tumors.  The  impact  of  epilepsy  in  the  US  is 
significant  with  a  total  cost  to  the  nation  for  seizures  and  epilepsy  of  approximately  $12.5  billion. 
Epilepsy  consists  of  more  than  40  clinical  syndromes  affecting  40  million  people  worldwide. 
Approximately  25  percent  of  individuals  receiving  antiepileptic  medication  have  inadequate 
seizure  control;  however,  80%  individuals  with  medication  resistant  epilepsy  might  be  cured 
through  surgery  if  one  were  able  to  precisely  localize  the  seizure  focus.  The  proposed  research 
will  significantly  advance  our  ability  to  localize  such  foci,  and  thereby  offer  curative  epilepsy 
surgery  for  this  devastating  disease.  Photoacoustic  tomography  (PAT)  uniquely  combines  the 
high  contrast  advantage  of  optical  imaging  and  the  high  resolution  advantage  of  ultrasound 
imaging  in  a  single  modality.  In  addition  to  high  resolution  structural  information,  the  proposed 
PAT  is  also  able  to  provide  functional  information  that  are  strongly  correlated  with  regional  or 
focal  seizure  activity,  including  blood  volume  and  blood  oxygenation  because  of  the  high 
sensitivity  of  optical  contrast  to  oxyhemoglobin  and  deoxyhemoglobin  concentrations.  The 
hypothesis  of  the  proposed  research  is  that  PAT  offers  the  possibility  to  non-invasively  track 
dynamical  changes  during  seizure  occurrence.  The  overall  goal  of  this  research  is  to  advance  a 
finite  element  based  photoacoustic  tomography  method  for  epilepsy  imaging,  using  both 
laboratory  and  in  vivo  experiments.  Specifically,  in  this  project  we  propose:  (1)  To  design, 
construct  and  test  a  transducer  array  system  for  both  2D  and  3D  PAT  imaging;  (2)  To  advance 
reconstruction  algorithms  and  associated  image  enhancement  schemes  for  quantitative  PAT;  (3) 
To  evaluate  and  optimize  the  integrated  functioning  of  the  hardware  and  software  components 
of  the  transducer  array-based  system,  using  simulation  and  phantom  experiments;  (4)  To  test 
and  validate  the  PAT  system  using  a  well  established  animal  model  of  temporal  lobe  epilepsy. 

Body 

This  report  describes  work  accomplished  in  Year  3  (Months  25-36)  of  the  project.  As  outlined  in 
the  approved  Statement  of  Work  (SOW),  the  tasks  during  this  period  of  time  include:  Task  5. 
Months  25-42:  Implement  reconstruction  codes  in  the  areas  of  advanced  boundary  conditions, 
3D  algorithms  and  parallelization  of  3D  codes;  perform  simulation  and  phantom  experiments  to 
evaluate  the  software  developments.  Task  6.  Months  25-42:  Conduct  animal  experiments  for 
evaluating  the  PAT  imaging  system. 

The  sections  below  consist  of  (1)  software  development,  and  (2)  animal  experiments  that  reflect 
the  tasks  associated  with  the  SOW  during  Months  25-36. 


1.  Software  Development  (Task  5) 

We  are  pleased  to  report  that  we  have  completed  the  software  development  proposed  for  Year 
3  including  advanced  boundary  conditions,  3D  algorithms,  parallelization  of  3D  codes,  and  their 
testing  using  simulation  and  phantom  experiments. 

1.1  Advanced  Boundary  Conditions 

1.1.1  Higher-order  absorption  boundary  conditions 
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Our  PAT  reconstruction  algorithm  uses  a  regularized  Newton  iterative  strategy  to  update  an 
initial  optical  or  acoustic  properties  distribution  to  minimize  an  object  function  composed  of  a 
weighted  sum  of  the  squared  difference  between  computed  and  measured  acoustic  data.  The 
computed  acoustic  data  is  obtained  by  solving  the  Helmholtz  photoacoustic  wave  equation 
using  finite  element  method.  In  the  frequency  domain,  the  photoacoustic  wave  equation 
(governing  equation)  can  be  written, 

V2p(r,  co)  +  kl  (1  +  0)p(r,  a>)  =  ik0  (1 ) 

where  p  is  the  pressure  wave;  k0  =colc 0  is  the  wave  number  described  by  the  angular 
frequency,  co  and  the  speed  of  acoustic  wave  in  a  reference  or  coupling  medium,  c0\  O  is  a 
coefficient  that  depends  on  both  acoustic  speed  and  attenuation. 

We  have  implemented  the  higher-order  absorption  boundary  conditions  to  improve  the  model 
accuracy  in  Eq.  (1).  We  also  compared  the  reconstructions  using  higher-order  and  low-order 
absorption  boundary  conditions.  For  low-order  absorption  boundary  conditions,  we  apply  the 
Sommerfeld  radiation  conditions  for  the  Helmholtz  wave  equation  on  the  boundary, 


-ik0p  =  Vp  ■  n  (2) 

For  the  higher  order  case,  a  second-order  absorption  boundary  conditions  for  arbitrarily  shaped 
convex  artificial  boundaries  of  3D  acoustic  problems  are  also  employed  here  for  comparison 
(Bamberger  1990), 


-Mp  =  Vp-n  on  27  (3) 

in  which  M  is  differential  operator  defined  as  the  function  of  the  mean  curvature  and  Gaussian 
curvature  at  a  point  of  a  surface,  27  is  the  triangulated  surfaces  of  a  3D  linear  tetrahedral 
element. 

To  test  the  influence  of  boundary  conditions  on  the  reconstruction  accuracy  and  algorithm 
performance,  3D  simulation  tests  are  conduced  here.  For  the  simulation  test  1,  one  6.0  mm- 
diameter  cylindrical  target  was  embedded  in  the  cylindrical  background  medium  (diameter: 
30.0mm;  height:  20.0mm),  as  shown  in  Fig.  1(a).  The  absorbed  light  energy  density  for  the 
background  and  target  were  0=  1  and  cp=3  (mJ/mm3),  respectively.  We  can  see  from  Figs.  1(b) 
and  1(c)  that  the  target  was  clearly  detected  by  using  both  the  low-order  and  higher-order 
absorption  boundary  conditions.  And  we  found  that  the  impact  of  the  boundary  conditions  on  the 
recovered  absorbed  energy  density  images  is  not  obvious  for  small  target  detection. 
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(c) 


Fig.  1  Reconstructed  absorbed  energy  density  images  for  the  one-target  case,  (a)  The  Test  geometry  for  the 
simulation  test;  (b)  both  3D  and  2D  slice  display  at  x=1.0mm  for  the  recovered  energy  density  with  higher-order 
absorption  boundary  conditions;  (c)  both  3D  and  2D  slice  display  at  x=1.0mm  for  the  recovered  energy  density  with 
low  order  absorption  boundary  conditions.  The  axes  (left  and  bottom)  illustrate  the  spatial  scale,  in  millimeters, 
whereas  the  color  scale  (right)  records  the  absorbed  optical  energy  density,  in  mJ/mm3. 

The  second  test  is  intended  to  test  the  influence  of  boundary  conditions  on  big  target  detection. 
In  this  case,  two  15  mm-diameter  cylindrical  objects  mimicking  two  seizure  focuses  were 
embedded  in  the  background  medium  where  the  spacing  between  the  two  “focuses”  was  2mm. 
The  absorbed  light  energy  density  for  the  background  and  “seizure  focuses”  were  0=  1  and  <P=3 
(mJ/mm3),  respectively.  We  can  see  from  Fig.  2  that  both  the  2  targets  can  be  clearly  detected 
by  using  both  the  low-order  and  high-  order  absorption  boundary  conditions.  However,  the 
impact  of  the  boundary  conditions  on  the  reconstruction  of  large  targets  is  obvious,  where 
strong  boundary  artifacts  are  observed  for  the  reconstruction  using  the  low-order  absorption 
boundary  conditions  (Fig.  2a). 
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Fig.  2  The  recovered  absorbed  energy  density  images  with  the  low-order  absorption  boundary  conditions  (a),  and  the 
higher-  order  absorption  boundary  conditions  (b). 

1.1.2  Mismatches  between  the  computational  model  and  the  physical  data  acquisition 

One  of  our  algorithmic  development  activities  will  be  on  issues  related  to  mismatches  between 
the  computational  model  and  the  physical  data  acquisition.  In  this  regard  we  have  investigated 
the  possible  model  representations  of  the  transducers  and  concomitant  boundary  conditions.  In 
particular,  we  implemented  a  model  that  considers  the  finite  ultrasound  beam  effects  associated 
with  the  receiver  hardware.  According  to  the  simulation  tests,  we  found  that  the  ultrasound 
beam  size  has  no  effect  on  the  final  reconstructed  results  when  our  finite  element  based 
reconstruction  method  is  used,  and  the  results  for  a  typical  test  are  provided  in  Fig.  3. 


(c) 

Fig.  3  Reconstructed  absorbed  energy  density  images  for  a  two-target  case,  (a)  The  test  geometry  for  the  simulation 
test;  (b)  both  3D  and  2D  slice  display  for  the  recovered  energy  density  using  point  detectors;  (c)  both  3D  and  2D  slice 
display  for  the  recovered  energy  density  using  finite  ultrasound  beam  sizes.  The  axes  (left  and  bottom)  illustrate  the 
spatial  scale,  in  millimeters,  whereas  the  color  scale  (right)  records  the  absorbed  optical  energy  density,  in  mJ/mm3. 
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1.2  3D  Algorithms  and  Parallelization  of  3D  Codes 
1.2.1  3D  Algorithms 

We  have  implemented  3D  algorithms  that  can  reconstruct  the  absolute  optical  absorption 
coefficient  distribution  of  3D  heterogeneous  medium  when  conventional  PAT  method  is 
combined  with  a  light  transport  model.  Therefore,  the  capabilities  of  PAT  have  now  been 
extended  to  provide  spatially  resolved  quantitative  physiological  and  molecular  information  that 
can  be  derived  from  tissue  absorption  spectra.  Basically,  the  reconstruction  method  includes 
two  steps.  The  first  is  to  obtain  the  3D  images  of  absorbed  light  energy  density  through  our  3D 
PAT  reconstruction  algorithm  that  is  based  on  the  finite  element  solution  to  the  photoacoustic 
wave  equation  in  frequency  domain  subject  to  the  radiation  or  absorbing  boundary  conditions. 
The  second  step  is  to  recover  the  optical  absorption  coefficient  distribution  using  the  photon 
diffusion  model  coupled  with  the  absorbed  optical  energy  density  images  obtained  from  the  first 
step. 

In  the  first  step,  our  3D  PA  reconstruction  algorithm  uses  a  regularized  Newton  iterative  strategy 
to  update  an  initial  absorbed  optical  energy  density  distribution  to  minimize  an  object  function 
composed  of  a  weighted  sum  of  the  squared  difference  between  computed  and  measured 
acoustic  data.  In  the  second  step,  to  recover  the  optical  absorption  coefficient  from  the 
absorbed  energy  density,  <6  ,  the  photon  diffusion  equation  as  well  as  the  Robin  boundary 
conditions  can  be  written  in  consideration  of  4>  =  , 


V  •  D(r)V  ( E(r)(S>(r ))  -  <&(r)  =  -S(r)  (4) 

-DV(£(r)0)-B  =  £(r)a®  (5) 

where  E(r)  =  l /Ma(r) .  D(r)  is  the  diffusion  coefficient  and  D  =  l/(3(jua+Ml))  where  /./'  is  the  reduced 
scattering  coefficients,  a  is  a  boundary  condition  coefficient  related  to  the  internal  reflection  at 
the  boundary,  and  S(r)  is  the  incident  point  or  distributed  source  term.  For  the  inverse 
computation,  the  so-called  Tikhonov-regularization  sets  up  a  weighted  term  as  well  as  a  penalty 
term  in  order  to  minimize  the  squared  differences  between  computed  and  measured  absorbed 
energy  density  values, 


min{||<Dc  -<P"||2  +  A||L[E-E0]||2  (6) 

where  L  is  the  regularization  matrix  or  filter  matrix,  f5  is  the  regularization  parameter, 
0>°  and  <DC  =  (®j,®c2,...,®^)T  where  ®"  is  the  absorbed  energy  density  obtained 

from  PAT,  and  ®j  is  the  absorbed  energy  density  computed  from  Eqs.  (4)  and  (5)  for  /=  1,  2..., 
N  locations  within  the  entire  PAT  reconstruction  domain.  The  initial  estimate  of  absorption 
coefficient  can  be  updated  based  on  iterative  Newton  method  as  follows  with  /?=1, 

A(E)  =  (JrJ  +  AI  +  LrL)_1[Jr (®  -  )]  (7) 

in  which  Laplacian-type  filter  matrix,  L  ,  is  employed  and  its  elements  are  written, 
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1 

if  i  =  j 

-1/AW 

if  i,j  a  one  region 

0 

if  i,j  cdiffrent  region 

(8) 


where  NN  is  the  total  node  number  within  one  region  or  tissue.  The  optical  absorption  coefficient 
distribution  is  then  reconstructed  through  the  iterative  procedures  described  by  Eqs.  (4)  and  (7) 
where  a  priori  spatial  information  from  the  PAT  image  is  incorporated  in  through  the  matrix  L. 

Phantom  experiments  were  conducted  to  test  the  3D  algorithms  using  the  192-channel  PAT 
system  developed  in  Years  1  and  2.  The  3D  phantom  was  made  of  Intralipid,  India  ink,  distilled 
water  and  Agar  powder.  For  the  phantom,  the  two  “seizure  focuses”  had:  (absorption  coefficient) 
Ha=  0.07mm'1,  (reduced  scattering  coefficient) /is  =2. 5mm'1  and  a  diameter  of  6mm,  and  the 

soft  tissues  between  the  seizure  focuses  had:  /u  =  0.015mm'1,  and  ju's  =1 .5mm"1. 


The  reconstructed  results  for  the  phantom  experiment  are  shown  in  Fig.  4,  where  Fig.  4a 
displays  the  coronal  cross  section  of  the  reconstructed  3D  absorbed  energy  density  image.  We 
can  observe  that  the  seizure  focuses  are  clearly  identified  (the  red  area).  The  recovered 
absorption  coefficient  image  is  displayed  in  Fig  4b,  where  the  targets  are  recovered  in  terms  of 
the  position,  size  and  quantitative  optical  property. 


(a)  (b) 

Fig.  4  (a)  The  selected  coronal  slice  from  the  recovered  absorbed  energy  density  image,  (b)  the  selected  coronal 
slice  from  the  recovered  absorption  coefficient  image. 


1.2.2  Parallelization  of  3D  reconstruction  codes 


As  we  have  experienced,  full  3D  reconstruction  requires  significant  computational  times 
especially  at  higher  ultrasound  frequencies.  Given  the  large  problem  size,  three  major  areas  of 
the  algorithm  are  costly:  (i)  computation  of  the  forward  solution,  (ii)  buildup  of  the  Jacobian 
matrix  which  requires  calculation  of  the  derivatives  of  acoustic  field  with  respect  to 
optical/acoustic  properties  at  the  receiver  sites,  and  (Hi)  computation  of  a  full  matrix  equation  for 
the  inverse  solution.  Because  of  the  massively  parallel  structure  of  these  types  of  image 
reconstruction  algorithms,  parallel  computing  was  adapted  for  our  3D  PAT  reconstruction  and 
the  code  was  developed  to  resolve  this  large-scale  computational  problem. 

Parallel  computing  is  basically  the  concurrent  execution  of  many  computations  using  many 
processors.  In  a  parallel  computing  scheme,  the  memory  request  and  the  computation  tasks  for 
a  computation  assignment  can  be  spread  on  multiple  processors  as  shown  in  Fig.  5.  In  the 
following,  we  will  give  a  detailed  description  on  the  parallel  computing  based  3D  PAT  algorithm. 
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In  the  finite  element  based  PAT  algorithm,  the  forward  calculation  of  acoustic  field  and  the 
determination  of  Jacobian  matrix  J  are  based  on  the  forward  matrix  A  at  each  frequency  , 

ordered  from  g\  to  coK  with  total  AT  frequency  elements.  The  matrix  A  in  forward  equations  is 
usually  a  symmetric  sparse  matrix  and  stored  in  memory  by  banded  storage  or  compressed 
storage  strategy,  requiring  less  memory  than  the  Jacobian  matrix  J  and  the  Hessian  matrix  JTJ  , 
which  is  usually  a  full  matrix.  Herein,  we  spread  the  forward  calculation  of  acoustic  field  and  the 
determination  of  Jacobian  matrix  on  distributed  processors  (total  number  is  Q  +  l)  by  the 

frequency  element ®.,  where  the  whole  matrix  A  under  certain  frequency  elements  is  stored  in 
the  memory  of  the  specifically  assigned  processor.  As  shown  in  flow  chart  of  the  high 
performance  photoacoustic  tomography  (Fig.  6),  the  matrix  A  at  frequency^,  is  stored  in  the 
memory  of  a  processor  determined  by  Mod{i-\,L )  ,  where  L  is  the  average  number  of 
frequency  elements  on  each  processor  (total  frequency  elements  K  divided  by  total 
processors  Q  +  \  ).  In  each  processor,  the  forward  modeled  acoustic  field  as  well  as  the 
elements  of  Jacobian  matrix  at  the  frequencies  associated  with  the  specified  processor  is 
calculated  independently,  after  which  the  Jacobian  matrix  J  over  the  whole  frequency  range  is 
assembled  and  stored  dispersedly  in  the  memories  of  the  distributed  processors.  In  this  way  of 
parallel  computing  and  storage,  the  computation  load  and  storage  is  evenly  assigned  among  the 
processors,  and  the  computation  assignment  is  maximally  parallelized  since  each  processor 
can  run  the  assigned  task  independently  without  communication  with  other  processors. 

Jacobian  matrix  J  and  the  Hessian  matrix  JTJ  (denoted  as  H )  are  full  matrices,  requiring 
distributed  storage  over  the  processors.  A  wrapping  storage  over  the  columns  of  the  Jacobian 
matrix  J  is  used  to  evenly  store  the  Jacobian  sub-matrix  over  the  processors  and  further 
efficiently  calculate  the  Hessian  matrix //with  minimally  mutual  communication  between  the 
processors.  The  wrapping  storage  is  described  below, 

{jj^eCPU  Myid,  if  Myid  =  Modij +  (9) 

where  =  ,j  =  l,...,N. 

With  the  wrapping  storage  of  the  Jacobian  matrix  J,  the  Hessian  matrix //can  be  calculated  and 
stored  by 

\HJ }  =  { J,  J2  ■  •  •  Jj } r  \jj  }  e  CPU  Myid,  if  Myid  =  Modij  - 1,  Q  + 1)  (10) 

where  {//.}  =  {z^.,...,/^.}' ,  j  =  . 

Because  the  Hessian  matrix  H  is  symmetric,  we  store  only  the  upper  triangle  of  the  Hessian 
matrix  spread  over  the  processors.  In  the  calculation  of  element//^,  communication  and  data 

exchange  may  be  involved  among  the  distributed  processors,  except  that  the  vector  jj;  j  and 

j  J;  |  are  both  stored  in  the  memory  of  the  same  processor.  Since  the  Hessian  matrix  is  stored 

dispersedly,  the  inverse  reconstruction  requires  a  parallel  solving  method,  where  parallelized 
Cholesky  decomposition  is  used  in  our  high  performance  PAT. 
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With  the  high  performance  PAT,  simulated  3D  blood  vessels  with  different  diameters  (0.6mm, 
0.28mm  and  0.14  mm)  can  be  accurately  reconstructed,  as  shown  in  Fig.  7  for  one  slice  of  the 
3D  image. 


A  Computation  Assignment 


l 


Fig-5  Schematic  of  parallel  computing  with  a  combination  of  both  shared  memory  and  distributed  memory. 
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Fig.  6  Flow  chart  of  the  high  performance  PAT  based  on  parallel  computing  technique  and  finite  element  method. 
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Fig.  7  The  reconstructed  coronal  slice  selected  from  simulated  3D  blood  vessels. 


2.  Animal  Experiments  (Task  6) 

In  Year  3  we  conducted  extensive  animal  experiments  to  evaluate  the  PAT  system 
constructed  in  Years  1  and  2.  We  were  able  to  image  the  functional  anatomy  of  a  focal 
seizure  noninvasively  and  in  real  time  at  a  spatial  resolution  of  ~150pm  (Figs.  8-11).  In 

this  real  time  PAT  system  (Fig.  8a),  light  from  a  Ti:  Sapphire  laser  tunable  from  690  to  950  nm 
was  delivered  to  the  cortical  brain  surface  through  an  optical  fiber.  Part  of  the  energy  of  each 
laser  pulse  was  detected  by  a  photodiode  for  calibration.  A  192-element  full-ring  transducer 
array  was  used  to  capture  the  photoacoustic  signals  generated  by  the  laser  light.  All  the  192 
channels  were  equipped  with  preamplifiers  and  secondary  stage  amplifiers  for  optimized  signal- 
to-noise  ratio.  A  3:1  electronic  multiplexer  coupled  with  64-channel  analog-to-digital  converters 
completed  the  192-channel  data  acquisition,  providing  an  imaging  speed  of  0.33  s/frame 
primarily  limited  by  the  10  Hz  laser  repetition.  The  spatial  resolution  of  PAT  is  inversely  related 
to  the  bandwidth  of  the  ultrasonic  transducer.  In  our  PAT  system  each  ultrasonic  detector  had  a 
5  MHz  central  frequency,  a  70%  nominal  bandwidth,  and  a  diameter  of  6  mm  (Blatek,  Inc.,  PA, 
USA),  resulting  in  a  spatial  resolution  of  150pm. 

We  first  tested  the  PAT  system  accuracy  by  imaging  a  living  rat  brain  with  both  the  scalp  and 
skull  intact.  Results  indicate  that  the  noninvasive  PAT  image  of  the  rat  cortical  vasculature 
correlated  with  the  anatomical  location  of  each  vessel  (Figs.  8b,  c).  The  middle  cerebral  artery, 
right  hemispheres,  left  hemispheres,  left  olfactory  bulbs,  and  right  olfactory  bulbs  were  clearly 
visualized  in  the  PAT  image.  Numbers  1-5  marked  in  the  image  indicate  that  the  micro-blood 
vessels  have  diameters  of  less  than  100  pm. 

We  then  tested  the  system  accuracy  in  identifying  a  focal  EEG  seizure  and  focal  EEG 
epileptogenic  spike  and  wave  discharges.  Adult  Sprague  Dawley  rats  were  experimentally 
prepared  to  demonstrate  focal  cortical  seizures.  10  pi  of  1.9  mM  of  the  GABA  antagonist 
bicuculline  methiodide  and  10  pi  of  saline  control  were  injected  into  the  left  and  right  parietal 
neocortex,  respectively.  The  region  of  interest  including  the  injection  site  and  remaining  skull 
surface  was  photoacoustically  scanned  continuously  at  time  beginning  immediately  before  the 
electrographic  seizure  onset  time,  during  the  seizure  time,  and  after  the  seizure  ended  (Fig.  9c). 
Whereas  a  significant  increase  of  absorption  is  seen  in  the  region  of  the  bicuculline  injection,  no 
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absorption  contrast  is  observed  in  the  region  of  the  saline  injection  or  surrounding  brain  areas 
(Fig.  9c).  These  results  suggest  that  the  increase  in  local  and  surrounding  brain  tissue 
absorption  was  mostly  due  to  the  local  bicuculline  induced  seizures  (Figs.  9d,  e).  These  results 
demonstrate  the  local  network  of  the  seizure  foci,  revealing  complex  interactions  among 
different  neuronal  groups  at  a  spatial  scale  of  ~  100  pm.  In  summary,  these  results  highlight  the 
accuracy  of  PAT  for  noninvasive  precise  localization  of  epileptic  foci  and  network. 

We  also  employed  PAT  to  see  if  we  could  identify  interictal  spike  and  wave  epileptiform 
discharges.  Interictal  discharges  are  often  seen  in  focal  epilepsy  and  represent  highly 
synchronized  events.  Spikes  often  correlate  with  the  spatial  context  of  an  epileptic  focus.  We 
recorded  PAT  imaging  from  a  region  of  ~2  *  3  mm  in  the  purple  area  (Marked  ROI  in  Fig.  10a). 
We  observed  a  clear  optical  absorption  change  directly  associated  with  de-oxygenation  at  a 
wavelength  of  755  nm  (Fig.  10c).  At  seizure  onset  at  Imin  6.267  s,  the  focus  measured  0.2053 
mm2  and  increased  to  0.5004  mm2  during  the  categorical  electrographic  seizure  spike  and  wave 
frequency  and  amplitude  increases  (Fig.  10b).  Corresponding  rate  changes  (0.33  s)  spike  and 
wave  discharges  and  PAT  images  were  observed  in  each  experimental  trial. 

To  understand  the  dynamic  interactions  within  the  seizure  functional  circuitry,  we  use  the 
frequency  domain  pairwise  Granger  analysis  to  analyze  the  PAT  time  series  data.  The  time 
series  data  was  selected  from  200  sets  of  PAT  images  sampled  at  3Hz  within  the  seizure  onset 
period,  and  at  each  time  point,  we  chose  nine  ROIs  (Fig.  Ilci)  and  used  the  averaged  optical 
absorption  within  each  region  for  the  network  analysis.  The  nine  sets  of  time  series  analysis 
were  classified  into  3  groups:  (1)  the  primary  ictal  onset  area  and  corresponding  4  surrounding 
regions  of  interest  (ROI);  (2)  the  primary  seizure  focus,  secondary  homotopic  foci  and  2  ROIs 
surrounding  the  primary  focus  in  the  contralateral  cortex;  and  (3)  the  middle  cerebral  artery, 
primary  seizure  focus,  and  secondary  homotopic  seizure  focus.  From  group  1  analysis,  we  can 
see  the  causal  influence  from  the  primary  focus  to  its  surrounding  ROIs  (Fig.  11b,  left).  We  also 
computed  the  zero-lag  correlation  between  the  primary  focus  and  its  4  surrounding  ROIs.  We 
found  that  between  the  primary  focus  and  ROI  2,  there  was  an  obvious  negative  correlation  (- 
0.68),  indicating  that  the  neural  activities  within  these  two  regions  changed  in  opposite  direction 
over  time.  Such  inverted  change  in  optical  signal  was  likely  due  to  the  opposite  oxygen  delivery 
or  blood  flow.  Figure  11c2  shows  the  Granger  analysis  for  group  2  data  where  we  see  that  the 
primary  and  mirror  foci  influenced  mutually  with  a  stronger  influence  of  primary  focus  to  the 
mirror  focus.  Moreover,  we  observed  that  the  primary  focus  could  also  influence  the  other  2 
ROIs  in  the  contralateral  hemisphere.  Finally,  from  group  3  analyses  (Fig.  11c3),  we  can  see 
that  the  primary  focus  strongly  influenced  the  hemodynamic  change  in  the  middle  cerebral 
artery,  which  in  turn  showed  influence  on  the  secondary  homotopic  foci. 
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Figure  8  Real-time  PAT  system  for  seizure  dynamics,  (a)  Schematic  of  the  real-time  PAT  system.  A 
192-element  full-ring  transducer  array  was  used  to  capture  the  PA  signal  during  seizure  onset,  (b) 
Noninvasive  PAT  imaging  of  a  rat  brain  in  vivo  with  the  skin  and  skull  intact.  MCA,  middle  cerebral  artery; 
RH,  right  hemispheres;  LH,  left  hemispheres;  LOB,  left  olfactory  bulbs;  ROB,  Right  olfactory  bulbs,  (c) 
Open-skull  photograph  of  the  rat  brain  surface  acquired  after  the  PAT  experiment.  Numbers  1-5  indicate 
the  corresponding  blood  vessels  in  the  PAT  image  and  rat  brain  photograph. 


12 


Huabei  Jiang,  PhD/Annual  Report 


Figure  9  Noninvasive  epileptic  foci  localization,  (a)  EEG  recordings  of  seizure  onset  after  bicuculline 
methiodide  injection,  (b)  Schematic  showing  the  location  of  BMI  injection  and  saline  injection  (control).  B, 
BMI  injection;  S,  saline  injection,  (c)  Reconstructed  PAT  image  before  (left)  and  after  (right)  the 
bicuculline  methiodide  injection.  The  BMI  and  saline  were  injected  into  the  left  and  right  parietal  neocortex, 
respectively.  A  significant  increase  of  optical  absorption  is  seen  in  the  region  of  the  bicuculline  methiodide 
injection,  while  no  absorption  contrast  is  observed  in  the  region  of  the  saline  injection  relative  to  its 
surroundings,  (d)  Series  of  PAT  images  from  three  representative  transverse  planes  parallel  to  the  skin 
surface  with  755  nm  wavelength.  Slices  5,  7,  and  9  are  the  images  obtained  3,  5  and  7  mm  under  the  skin, 
respectively.  The  images  show  different  spatial  patterns  at  the  seizure  foci  at  different  tomographic  layers, 
(e)  Three  dimensional  rendering  of  the  epileptic  foci  from  different  tomographic  layers. 
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Figure  10  Real-time  monitoring  of  ictal  onset,  (a)  A  PAT  image  showing  the  cerebral  vasculature  and 
the  position  of  BMI  injection.  ROI,  region  of  interest;  Scale  bar:  2mm.  (b)  EEG  recordings  showing  the 
seizure  onset  about  1  min  after  BMI  injection,  (c)  PAT  images  recording  the  ictal  onset  in  real  time.  At 
1  min  6.267  s  the  area  of  the  focus  was  0.2053  mm2.  The  area  of  the  focus  then  increased  in  size  over  the 
next  2  s  to  0.5004mm2,  corresponding  to  an  increase  in  the  amplitude  of  the  EEG  spikes.  The  area  of  the 
seizure  focus  was  derived  from  the  PAT  images  by  thresholding  to  a  pixel  value  one  standard  deviation 
above  the  pixel  values  from  the  area  of  the  focus  during  the  control  conditions.  Scale  bar:  500pm. 
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Figure  11  In  vivo  Maps  of  contralateral  homotopic  seizure  foci  (a)  Left,  the  PAT  image  of  the  rat  brain; 
the  yellow  dotted  rectangular  shows  the  region  of  interest  for  analysis;  Right,  noninvasive  imaging  of  the 
primary  and  mirror  foci  in  contralateral  homotopic  cortex.  A1-4  shows  how  secondary  epileptic  foci  were 
generated  and  diminished.  The  PA  signal  in  the  secondary  foci  was  smaller  in  magnitude  and  delayed  in 
time  compared  to  the  signal  recorded  from  the  primary  foci,  (b)  Left,  Granger  analysis  of  the  image 
around  the  foci  and  its  immediate  secondary  foci;  Right,  PA  signal  from  the  foci  (red)  and  secondary  area 
(blue)  during  an  ictal  event,  (c)  Pairwise  Granger  analysis  of  the  seizure  foci,  contralateral  homotopic  foci 
and  middle  cerebral  artery.  Cl,  PAT  image  showing  the  9  selected  regions  of  interest;  C2,  Arrows 
indicate  Granger  connectivity  maps  between  the  primary  and  secondary  foci;  C3,  Connectivity  maps 
between  the  primary  and  secondary  foci,  and  the  middle  cerebral  artery. 

We  also  found  a  spatial  temporal  correlation  between  the  occurrence  of  interictal  spikes  and 
local  blood  vessel  dynamic  changes  of  the  vessel  size.  Figure  12a  shows  PA  image  of  the  rat 
cortex  vasculature  through  the  intact  scalp  and  skull.  The  red  arrow  indicates  the 
microvascularity  (MV)  along  the  direction  marked  by  the  yellow  dashed  line,  representing  a 
cross-sectional  scanning  using  a  50  MHz  ultrasound  transducer.  Positive  and  negative  acoustic 
peaks  induced  by  a  25  pm  blood  vessel  were  sorted  out  as  target  signals  to  track  the  change  of 
the  vessel  diameter.  Typical  PA  signals  from  the  targeted  vessel  at  different  times  show 
apparent  vessel  vasomotion  by  -40%  during  the  interictal  discharges  (Fig.  12b).  Interestingly, 
the  EEG  recording  (Fig.  12c)  shows  that  the  interictal  spikes  had  a  strong  correlation  with  the 
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discrete  change  in  vessel  size  observed  photoacoustically  (Fig.  12d;  error  bars  ±  s.d.  were 
calculated  from  20  consecutive  spikes).  Vasomotion  during  seizure  onset  and  spread  was 
clearly  caused  by  the  oscillation  in  vessel  diameter  and  the  expansion  of  the  vessel  cross 
section.  We  reasoned  that  the  dilation  resulted  in  an  active  cerebral  microvasculature  response 
to  increased  metabolic  activity  accompanying  the  seizure.  The  combination  of  multi-scale 
imaging  ability  and  endogenous  hemoglobin  contrast  makes  PAT  a  promising  tool  for  imaging 
microvasculature  during  seizure  onset. 


Figure  12  Changes  of  Blood  Vessel  Diameter  during  interictal  onset,  (a)  PA  image  of  the  rat  cortex 
vasculature  with  the  intact  scalp  and  skull.  The  red  arrow  indicates  the  microvascularity  (MV)  along  the 
yellow  dashed  line  direction  for  a  cross-sectional  scanning  using  a  50  MHz  ultrasound  transducer,  (b) 
Typical  photoacoustic  signals  from  the  targeted  vessel  at  different  times,  (c)  EEG  recordings  show  the 
interictal  spikes,  (d)  Change  of  the  vessel  size  captured  by  PAT.  There  was  a  clear  correlation  between 
the  interictal  spikes  and  the  changes  of  the  vessel  size.  Error  bars  (±  s.d.)  were  calculated  from  20 
consecutive  spikes. 


Key  Research  Accomplishments 
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1.  We  have  developed  the  advanced  boundary  conditions  scheme,  and  implemented  the  3D 
reconstruction  codes  and  their  parallelization  which  allowed  efficient  high  performance  PAT 
image  reconstruction. 

2.  We  have  conducted  simulation  and  phantom  experiments  that  confirmed  our  software 
enhancement. 

3.  We  have  performed  extensive  in  vivo  experiments  using  the  rat  epilepsy  model  to  evaluate 
the  PAT  system  we  developed  in  this  project.  The  results  show  that  we  were  able  to  image 
the  functional  anatomy  of  epileptic  foci  noninvasively  and  in  real  time  at  a  spatial 
resolution  of  ~1 50pm. 

Reportable  Outcomes  (see  the  Appendix  to  this  Summary  Report) 

We  have  published  two  journal  papers  during  Year  3.  In  addition,  we  expect  that  the  results 
generated  from  the  simulation,  phantom  and  animal  experiments  completed  in  Year  3  will  allow 
us  to  produce  numerous  journal  publications  and  conference  proceedings  in  Year  4. 

Conclusions 

We  have  made  a  significant  progress  that  has  fulfilled  the  statement  of  work  proposed  for  Year 
3  of  this  project.  The  real-time  PAT  system  completed  in  Years  1  and  2  has  provided  us  a 
platform  for  performing  extensive  animal  experiments.  In  Year  4  we  will  focus  on  animal  studies 
of  continuous  recording  of  EEG  and  PAT  over  a  long  period  of  time  to  monitor  seizure  activity 
continuously  to  characterize  changes  during  interictal,  ictal  and  post-ictal  periods  when  a 
spontaneous  seizure  occurs.  We  are  confident  that  we  will  be  able  to  fulfill  or  exceed  the  work 
statement  for  Year  4. 

Appendix 

(1)  L.  Yao,  H.  Jiang,  “Enhanced  photoacoustic  tomography  using  total  variation  minimization”, 
Applied  Optics  50,  5031-5041  (201 1 ). 
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A  total  variation  minimization  (TVM)-based  finite  element  reconstruction  algorithm  for  photoacoustic 
(PA)  tomography  is  described  in  this  paper.  This  algorithm  is  used  to  enhance  the  quality  of  reconstructed 
PA  images  with  time-domain  data.  Simulations  are  conducted  where  different  contrast  levels  between 
the  target  and  the  background,  different  noise  levels,  and  different  sizes  and  shapes  of  the  target  are 
studied  in  a  30  mm  diameter  circular  heterogeneous  background.  These  simulated  results  show  that 
the  quality  of  the  reconstructed  images  can  be  improved  significantly  due  to  the  decreased  sensitivity 
to  noise  effect  when  the  TVM  is  included  in  the  reconstruction  algorithm.  The  enhancement  is  further 
confirmed  using  experimental  data  obtained  from  several  phantom  experiments  and  an  in  vivo  animal 
experiment.  ©  2011  Optical  Society  of  America 
OCIS  codes:  100.2980,  170.3010,  170.5120,  170.6960. 


1.  Introduction 

Photoacoustic  tomography  (PAT)  is  an  emerging 
noninvasive  imaging  technique  that  combines  the 
merits  of  high  optical  contrast  and  high  ultrasound 
resolution  in  a  single  modality  [1—3].  In  PAT,  the 
Helmholtz-like  photoacoustic  (PA)  wave  equation 
has  been  commonly  used  as  an  accurate  model  for  de¬ 
scribing  laser-induced  acoustic  wave  propagation  in 
tissue.  While  analytical  reconstruction  methods  have 
been  used  for  PA  image  reconstructions,  the  finite 
element  method  (FEM)-based  approach  appears  to 
be  particularly  powerful  in  this  regard  [4,5] .  The  ad¬ 
vantages  of  the  FEM-based  PAT  method  include 
(1)  quantitative  imaging  capability  by  recovering  op¬ 
tical  absorption  coefficient,  (2)  elimination  of  the  as¬ 
sumption  of  homogeneous  acoustic  medium  needed 
in  analytical  methods,  (3)  accommodation  of  object 
boundary  irregularity,  and  (4)  appropriate  boundary 
conditions  implementations. 

Our  finite  element  PAT  approaches  are  based  on  a 
regularized  iterative  Newton  method,  which  have 
been  tested  and  evaluated  using  extensive  simula- 
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tions  and  phantom  experiments  in  both  the  fre¬ 
quency  and  time  domain  [4—7] .  Compared  with  the 
frequency  domain  method,  the  time  domain  method 
can  reconstruct  images  with  less  background  arti¬ 
facts,  especially  when  the  target  size  is  very  small, 
and  can  provide  a  more  accurately  recovered  target 
size  [7].  While  these  results  are  encouraging  and  pro¬ 
mising,  we  realize  that  measurement  noises  are  still 
the  major  factor  affecting  the  quantitative  accuracy 
of  the  reconstructed  images.  For  example,  errors  for 
quantitatively  recovering  optical  absorption  coeffi¬ 
cient  can  be  as  large  as  10%  for  simulated  data  with 
5%  noise  added  and  20%  for  experimental  data  [6,7] . 

In  an  effort  to  reduce  the  noise  effect  and  further 
enhance  the  quantitative  accuracy  of  PA  image 
reconstruction,  in  this  paper,  we  consider  a  total 
variation  minimization  (TVM)  scheme  within  our 
FEM-based  reconstruction  framework.  Our  existing 
reconstruction  algorithms  are  based  on  the  least- 
squares  criteria  (i.e.,  the  regularized  Newton  meth¬ 
od)  [8,9]  that  stand  on  the  statistical  argument  that 
the  least-squares  estimation  is  the  best  estimator 
over  an  entire  ensemble  of  all  possible  pictures.  Total 
variation,  on  the  other  hand,  measures  the  oscilla¬ 
tions  of  a  given  function  and  does  not  unduly  punish 
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discontinuities  [10,11],  Hence,  one  can  hypothesize 
that  a  hybrid  of  these  two  minimization  schemes 
should  be  able  to  provide  higher-quality  image  re¬ 
constructions.  In  fact,  the  strategy  of  finding  minimal 
total  variation  has  proved  to  be  successful  in  applica¬ 
tions  including  electrical-impedance  tomography 
[11],  microwave  imaging  [12],  image  processing 
[10,13,14],  optimal  design  [15],  and  diffuse  optical  to¬ 
mography  [16],  The  TVM  has  also  been  implemented 
and  tested  in  a  non-FEM-based  PA  reconstruction 
framework  [17-19], 

In  Section  2,  we  describe  in  detail  the  implementa¬ 
tion  of  TVM  within  our  existing  PAT  reconstruction 
framework  in  the  time  domain.  In  Section  3,  we  de¬ 
monstrate  the  enhanced  reconstruction  algorithm 


using  simulated,  phantom,  and  in  vivo  data.  The 
conclusions  are  made  in  Section  4 

2.  TVM  Scheme 

To  describe  our  TVM  method,  we  first  briefly  intro¬ 
duce  the  FEM-based  regularized  Newton  method  in 
the  time  domain.  The  time  domain  PA  wave  equation 
in  tissue  can  be  described  as  follows  [20]: 

„  ^  1  d2p(r,t)  <>(r )fidJ(t) 

=  -  Cp  dt  •  (1) 

where  p  is  the  pressure  wave,  u0  is  the  speed  of  acous¬ 
tic  wave  in  the  medium,  [S  is  the  thermal  expansion 
coefficient,  Cp  is  the  specific  heat,  O  is  the  absorbed 
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Fig.  1.  (Color  online)  Reconstructed  absorbed  energy  density  images  from  simulated  data  with  and  without  the  TVM  enhancement  under 
different  noise  levels  (case  1).  (a)  Without  the  TVM,  0%  noise,  (b)  with  the  TVM,  0%  noise,  (c)  without  the  TVM,  10%  noise,  (d)  with  the 
TVM,  10%  noise,  (e)  without  the  TVM,  25%  noise,  (f)  with  the  TVM,  25%  noise.  The  axes  (left  and  bottom)  illustrate  the  spatial  scale  in 
millimeters,  whereas  the  color  scale  (right)  records  the  absorbed  energy  density  in  millijoules  per  cubed  millimeter. 
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energy  density,  and  J(t)=8(t-t0)  is  assumed  in 
our  study 

Expanding  p  as  the  sum  of  coefficients  multiplied 
by  a  set  of  basis  function  y/j,  p  =  ]T  pji//j,  the  finite 


Fig.  2.  (Color  online)  Comparison  of  the  exact  and  reconstructed 
absorbed  energy  density  profiles  along  transect  y  =  0  mm  for  the 
images  appearing  in  Fig.  1.  (a)  0%  noise,  (b)  10%  noise,  (c)  25% 
noise. 


element  discretization  of  Eq.  (1)  can  then  be  written 
as  [7] 


The  first-order  absorbing  boundary  conditions 
used  are  [21] 


Vp  •  h 


1  dp 
v0  d t 


P_ 
2  r’ 


(3) 


where  h  is  the  unit  normal  vector. 

In  both  the  forward  and  inverse  calculations,  the 
unknown  coefficient  <t>  needs  to  be  expanded  in  a 
similar  fashion  top  as  a  sum  of  unknown  parameters 
multiplied  by  a  set  of  locally  spatially  varying 
Lagrange  polynomial  basis  functions.  Thus,  the  ma¬ 
trix  form  of  Eq.  (2)  becomes 


[K}{p}  +  [C]{p}  +  [M]{p}  =  {B},  (4) 


where  the  elements  of  matrix  [. K\ ,  [C],  and  [M]  are 


and  the  column  vectors  {p},  {p},  {p},  and  {B\  are 


A  =  ^/s«(p''A)dS'w 

{p}  =  {piiP2,---pn}t; 

{p}  =  {PiiP2,---Pn}T-, 

{p}  =  {PiiP2,---Pn}T- 


Here,  the  Newmark’s  time-stepping  scheme  has  been 
used  for  the  discretization  of  time  dimension  [22,23], 
which  is  a  commonly  used  implicit  method  for  the 
second-order  propagation  equations  such  as  Eq.  (4). 

To  form  an  image  from  a  presumably  uniform  in¬ 
itial  guess  of  the  absorbed  energy  density  distribu¬ 
tion,  we  need  a  method  of  updating  O  from  its 
starting  value.  This  update  is  accomplished  through 
the  least-squares  minimization  of  the  following 
functional: 


M 

E(p,0)  =  ^(p°-pp2,  (5) 

j= i 

where  pj  and  pj  are  observed  and  computed  acoustic 
field  data  for  j  =  1, 2,  •  •  ■  M  boundary  locations.  Using 
the  regularized  Newton  method,  we  obtained  the 
following  equation  for  updating  O: 
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Fig.  3.  (Color  online)  Reconstructed  absorbed  energy  density  images  from  simulated  data  with  and  without  the  TVM  enhancement  with 
different  target  sizes  (case  2).  (a)  Without  the  TVM,  2  mm  diameter  target;  (b)  with  the  TVM,  2  mm  diameter  target;  (c)  absorbed  energy 
density  profiles  along  the  transect  y  =  0  mm  for  the  images  appearing  in  Figs.  1(a)  and  1(b)  (4  mm  diameter  target);  (d)  absorbed  energy 
density  property  profiles  along  the  transect  y  =  0  mm  for  the  images  appearing  in  Figs.  3(a)  and  3(b)  (2  mm  diameter  target).  In  (a)  and  (b), 
the  axes  (left  and  bottom)  illustrate  the  spatial  scale  in  millimeters,  whereas  the  color  scale  (right)  records  the  absorbed  energy  density  in 
millijoules  per  cubed  millimeter. 


(S7,S  +  2I)A/  =  ST(p0-^),  (6) 

where  p°  =  (p\,p°2,  ■  ■  -phY,  pc  =  (pcvp2 ,  •  •  ■pcM)T,  Ay 
is  the  update  vector  for  the  absorbed  optical  energy 
density,  S'  is  the  Jacobian  matrix  formed  by  dp/dO  at 
the  boundary  measurement  sites,  A  is  the  regulariza¬ 
tion  parameter  determined  by  combined  Marquardt 
and  Tikhonov  regularization  schemes,  and  I  is  the 
identity  matrix. 

Two  typical  approaches  exist  for  minimizing  total 
variation:  a  constrained  minimization  through  the 
solution  of  the  nonlinear  PA  equation  [8,24]  and  an 
unconstrained  minimization  by  addition  of  the  total 
variation  as  a  penalty  term  to  the  least-squares  func¬ 
tional  [10,13,14,25].  From  a  computational  stand¬ 
point,  unconstrained  minimizations  are  much 
easier  to  implement  and  require  fewer  modifications 
to  the  existing  algorithm  [13],  Thus,  in  this  study,  the 
unconstrained  TVM  is  used. 

We  now  incorporate  the  total  variation  of  O  as  a 
penalty  term  by  defining  a  new  functional  [9,10,13]: 

F(p,0)=F(p,<P)+L(0).  (7) 

Here, 


L(0)=  /  J<4|VO|2+«52dxdy 


(8) 


is  the  penalty  term,  and  cu<j>  and  <5  are  typically 
positive  parameters  that  need  to  be  determined  nu¬ 
merically.  The  minimization  of  Eq.  (7)  proceeds  in 
standard  fashion  by  the  differentiation  of  F  with  re¬ 
spect  to  each  nodal  parameter  that  constitutes  the  O 
distribution;  simultaneously,  all  these  relations  are 
set  to  zero.  These  steps  lead  to  the  following  non¬ 
linear  system  of  equations 


3 F_ 
40, 


M 


j=  1 


dpi 

®W,  +  Vt’ 


(i  =  1,2  •  •  -N), 


where 


(9) 


x  (body. 
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Fig.  4.  (Color  online)  Reconstructed  absorbed  energy  density  images  from  simulated  data  with  and  without  the  TVM  enhancement  with 
different  contrast  levels  between  the  target  and  the  background  (case  3).  (a)  Without  the  TVM,  1.5 1 1  contrast;  (b)  with  the  TVM,  1.5 : 1 
contrast;  (c)  absorbed  energy  density  profiles  along  the  transect  y  =  0  mm  for  the  images  appearing  in  Figs.  1(a)  and  1(b)  (2 : 1  contrast); 
(d)  absorbed  energy  density  profiles  along  the  transecty  =  0  mm  for  the  images  appearing  in  Figs.  4(a)  and  4(b)  (1.5 : 1  contrast).  In  Figs.  4 
(a)  and  4(b),  the  axes  (left  and  bottom)  illustrate  the  spatial  scale  in  millimeters,  whereas  the  color  scale  (right)  records  the  absorbed  energy 
density  in  millijoules  per  cubed  millimeter. 


Similar  to  Eq.  (6),  the  following  matrix  equation  for 
TVM  constrained  inversion  can  be  obtained: 

(STS  +  R+AI)AX  =  ST(p°-pc)-V,  (10) 

where 


avj  avj 

avj 

a®!  a®2 

'a®w 

av2  av2 .  _ 

av2 

a®!  a®2 

a®w 

dVNdVN _ 

dVN 

a®!  a®2 

and  V={V1,V2,---Vn}t. 

3.  Results  and  Discussion 

In  this  section,  the  TVM-enhanced  reconstruction 
algorithm  is  tested  and  evaluated  using  both 
simulated  and  experimental  data.  For  comparative 
purposes,  reconstructions  without  the  TVM  enhance¬ 
ment  are  also  presented. 


A.  Simulations 

In  these  simulations,  a  dual-meshing  method,  as  de¬ 
tailed  elsewhere  [6],  is  used  for  fast  yet  accurate  in¬ 
verse  computation.  The  fine  mesh  used  for  the 
forward  calculation  consisted  of  3627  nodes  and 
7072  elements,  while  the  coarse  mesh  used  for  the 
inverse  calculation  had  930  nodes  and  1768  ele¬ 
ments.  All  the  images  obtained  from  the  method 
without  the  TVM  are  the  results  of  three  iterations, 
while  those  obtained  from  the  TVM-enhanced  meth¬ 
od  are  the  results  of  15  or  more  iterations.  Parallel 
code  was  used  to  perform  these  calculations  on 
Beowulf  clusters  with  8  central  processing  units 
(CPUs).  The  computational  time  can  be  further  re¬ 
duced  if  clusters  with  more  than  8  CPUs  are  used.  As 
mentioned  earlier,  the  parameters  a>$,  and  <5  were 
determined  through  numerical  experimentation.  A 
constant  value  of  8  =  0.001  was  sufficient  for  the  cur¬ 
rent  simulation  and  experimental  studies,  while  the 
value  of  «,[>  is  related  to  the  signal-to-noise  ratio  of 
the  measurements  [14],  For  the  simulations  pre¬ 
sented,  a»<D  =  0.5  for  cases  1,  2,  and  3,  and  co $  = 
1.0  for  case  4  was  used. 

For  the  first  simulation,  the  test  geometry  is  a 
30  mm  diameter  circular  background  containing  an 
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Fig.  5.  (Color  online)  Reconstructed  absorbed  energy  density 
images  from  simulated  data  with  and  without  the  TVM  enhance¬ 
ment  for  three  targets  having  different  shapes  (case  4).  (a)  Exact 
image,  (b)  without  the  TVM,  (c)  with  the  TVM,  (d)  absorbed  energy 
density  profiles  along  the  transect  y  =  0  mm.  In  (a)-(c),  the  axes 
(left  and  bottom)  illustrate  the  spatial  scale  in  millimeters, 
whereas  the  color  scale  (right)  records  the  absorbed  energy  density 
in  millijoules  per  cubed  millimeter. 


off-center  4  mm  diameter  target  region.  The  target 
had  O  =  2.0  mJ/mm3,  while  the  background  had 
<t>  =  1.0  mJ/mm3.  In  this  case,  the  image  reconstruc¬ 
tion  was  performed  with  0%,  10%,  and  25%  noise 
added  to  the  “measured”  data.  Figure  1  gives  two  sets 
of  absorbed  energy  density  images  recovered  using 
the  reconstruction  method  without  [Figs.  1(a),  1(c), 
and  1(e)]  and  with  the  TVM  enhancement  [Figs.  1(b), 
1(d),  and  1(f)]  under  the  conditions  of  different  noise  le¬ 
vels.  As  can  be  seen,  enhancement  of  the  reconstruc¬ 
tion  by  the  incorporation  of  the  TVM  is  obvious  over 
the  method  without  the  TVM.  To  provide  a  more  quan¬ 
titative  assessment  of  these  images,  Fig.  2  is  included, 
in  which  the  reconstructed  absorbed  energy  density 
profiles  are  displayed  along  transects  through  the 
target  for  the  images  shown  in  Fig.  1.  We  find  that 
the  TVM-enhanced  method  not  only  can  improve 
the  quality  of  the  recovered  images,  but  also  can  de¬ 
crease  the  sensitivity  of  the  method  to  noise  effect. 

The  second  simulation  is  intended  to  investigate 
how  the  target  size  affects  the  TVM  enhancement. 
In  this  case,  no  noise  was  added  to  the  “measured” 
data,  and  the  diameter  of  the  off-center  target  was 
2  mm.  Figures  3(a)  and  3(b)  give  two  sets  of  absorbed 
energy  density  images  recovered  using  the  method 
without  the  TVM  [Fig.  3(a)]  and  with  the  enhance¬ 
ment  [Fig.  3(b)],  while  Figs.  3(c)  and  3(d)  present 
the  quantitative  profiles  of  the  absorbed  energy  den¬ 
sity  along  transects  through  the  target  for  the  images 
shown  in  Figs.  1(a),  1(b),  3(a),  and  3(b).  Again,  con¬ 
siderable  improvement  can  be  observed  from  the 
reconstructed  results  when  the  TVM  is  invoked  com¬ 
pared  with  the  method  without  the  TVM.  It  is  also 
interesting  to  note  that  the  enhancement  for  this 
case  is  more  striking  than  that  for  the  first  one  where 
the  background  contained  a  larger  target. 

The  third  simulation  aims  to  see  how  the  contrast 
level  ofthe  absorbed  energy  density  between  the  target 
and  the  background  impacts  the  TVM  enhancement. 
In  this  case,  no  noise  was  added  to  the  “measured” 
data,  and  the  off-center  target  had  a  diameter  of 
4  mm  and  O  =  1.5  mJ/mm3.  Figures  4(a)  and  4(b)  pre¬ 
sent  two  sets  of  absorbed  energy  density  images  recov¬ 
ered  using  the  method  without  the  TVM  [Fig.  4(a)]  and 
with  the  TVM  enhancement  [Fig.  4(b)]  >  while  Figs.  4(c) 
and  4(d)  show  the  quantitative  profiles  of  the  absorbed 
energy  density  along  transects  through  the  target  for 
the  images  shown  in  Figs.  1(a),  1(b),  4(a),  and  4(b).  We 
can  see  that  the  images  formed  by  incorporation  of  the 
TVM  are  clearly  enhanced  qualitatively  in  visual  con¬ 
tent  relative  to  that  obtained  using  the  method  with¬ 
out  the  TVM.  We  also  note  that  lower  the  contrast  level 
is,  the  larger  the  enhancement. 

In  the  fourth  simulation,  three  targets  having  dif¬ 
ferent  shapes  (circular,  elliptical,  and  rectangular) 
were  embedded  in  the  background,  and  the  absorbed 
energy  density  of  these  targets  was  2.5,  1.5,  and 
2.0 mJ/mm3,  respectively.  A  noise  level  of  25%  was 
added  to  the  “measured”  data  in  this  case.  Figure  5 
shows  the  exact  and  the  recovered  absorbed  energy 
density  images  as  well  as  the  quantitative  absorbed 
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Fig.  6.  UQI  calculated  from  the  recovered  images  with  and  without  TVM  enhancement  from  simulated  data,  (a)  Case  1  when  different 
noise  levels  are  considered,  (b)  cases  2-4. 


Fig.  7.  (Color  online)  Reconstructed  absorbed  energy  density  images  from  the  three  phantom  experiments,  (a)  Case  1  without  the 
TVM,  (b)  case  1  with  the  TVM,  (c)  case  2  without  the  TVM,  (d)  case  2  with  the  TVM,  (e)  case  3  without  the  TVM,  (f)  case  3  with  the 
TVM.  The  axes  (left  and  bottom)  illustrate  the  spatial  scale  in  millimeters,  whereas  the  color  scale  (right)  records  the  absorbed  energy 
density  in  millijoules  per  cubed  millimeter. 
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energy  density  property  profiles  along  the  transect 
that  across  these  targets.  Again,  the  improvement 
in  image  quality  is  apparent. 

In  order  to  quantitatively  evaluate  the  reconstruc¬ 
tion  quality  using  the  reconstruction  method  with 
and  without  TVM  enhancement,  we  used  the  univer¬ 
sal  quality  index  (UQI)  to  measure  the  degree  of  si¬ 
milarity  between  the  reconstructed  and  exact  images 
[26].  We  first  interpreted  an  image  f  as  vectors  of  size 
N  :  f  =  (f 21  • ' '  f n)T i  where  superscript  T  denotes 
the  matrix  transpose,  and  N  denotes  the  number 
of  image  data  acquired  from  the  FEM-based  algo¬ 
rithm.  We  then  defined  image  means,  variances  and 
covariances  over  the  whole  image  domain  as 

an 

V  k=l 


<jj  = 


B/t-f-o2, 


k=l 


(12) 


where  j  =  0  and  1,  and 


Covjf1^0}  =  ~f1)(fk  ~f°)- 


(13) 


k=i 


Hence,  the  UQI  is  defined  as 

2Cov{f1,f0}  2f1f° 


UQI^f0} 


(^)2+K)2r)2+(/°) 


f0\2  ' 


(14) 


UQI  measures  the  image  similarity  between  the  re¬ 
constructed  (f1)  and  reference/exact  (f°)  images,  and 
its  value  ranges  between  0  and  1.  The  value  of  the 
UQI  is  closer  to  1  when  the  reconstructed  image  is 
more  similar  to  the  exact  image.  We  calculated  the 
UQIs  for  the  simulation  cases  presented  earlier, 
and  the  results  are  given  in  Fig.  6,  where  Figs.  6(a) 
and  6(b)  show  the  UQIs  for  case  1  when  different  noise 
level  is  considered  and  for  all  the  four  cases  without 
added  noise.  In  Fig.  6,  the  horizon  axis  indicates  the 
case  number  (from  cases  1  to  4,  presented  earlier),  and 


X  position  (mm) 


Fig.  8.  (Color  online)  Recovered  absorbed  energy  density  profiles  along  (a)  y  =  -7.0  mm  crossing  the  3  mm  diameter  target  for 
experimental  case  1,  (b)  y  =  8.0  mm  crossing  the  2  mm  diameter  target  for  experimental  case  1,  (c)  y  =  1.0  mm  for  experimental  case 
2,  and  (d)  y  =  6.5  mm  for  experimental  case  3. 
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Fig.  9.  (Color  online)  CNR  calculated  for  the  recovered  images 
using  the  method  with  and  without  TVM.  (a)-(c)  Images  shown 
in  Figs.  7(b),  7(d),  and  7(f)  with  the  ROIs  marked  (four  pairs  of 
ROIs:  solid  line  circle,  t-ROI;  dashed  line  circle,  6-ROI).  (d)  CNR 
computed  for  the  four  pairs  of  ROIs  shown  in  (a)-(c).  In  Figs.  9 
(a)-9(c),  the  axes  (left  and  bottom)  illustrate  the  spatial  scale  in 
millimeters,  whereas  the  color  scale  (right)  records  the  absorbed 
energy  density  in  millijoules  per  cubed  millimeter. 


the  vertical  axis  shows  the  value  of  UQI.  We  immedi¬ 
ately  note  from  Fig.  6  that  the  TVM-based  method 
provides  significantly  better  UQIs  than  that  without 
TVM. 

B.  Phantom  and  In  Vivo  Experiments 

In  this  section,  both  phantom  and  in  vivo  experimen¬ 
tal  data  were  used  to  confirm  our  findings  from  the 
simulations.  The  experimental  setup  used  for  collect¬ 
ing  the  phantom  data  was  a  pulsed  Nd:YAG  laser- 
based  single  transducer  (1MHz)  scanning  system, 
which  was  described  in  detail  elsewhere  [5],  Three 
phantom  experiments  were  conducted.  In  the  first 
two  experiments,  we  embedded  one  or  two  objects 
with  a  size  ranging  from  3  to  0.5  mm  in  a  50  mm  dia¬ 
meter  solid  cylindrical  phantom.  The  phantom  mate¬ 
rials  used  consisted  of  Intralipid  as  a  scatterer  and 
India  ink  as  an  absorber  with  Agar  powder  (1%— 
2%)  for  solidifying  the  Intralipid  and  India  ink  solu¬ 
tion.  The  absorption  coefficient  of  the  background 
phantom  was  0.01mm-1,  while  the  absorption 
coefficient  of  the  target(s)  was  0.03  mm-1.  In  the  last 
experiment,  we  used  a  single  target-containing 
phantom,  aiming  to  test  the  capability  of  detecting 
targets  having  low  optical  contrasts  relative  to  the 
background  phantom.  In  this  case,  the  target  had 
an  absorption  coefficient  of  0.015  mm-1.  The  reduced 
scattering  coefficients  of  the  background  phantom 
and  targets  were  1.0  and  3.0  mm-1  for  the  first  two 
experiments,  and  1.0  and  2.0  mm-1  for  the  last 
experiment. 

A  total  of  120  receivers  were  equally  distributed 
along  the  surface  of  the  circular  background  region. 
In  the  reconstructions,  the  fine  mesh  used  for  the  for¬ 
ward  calculation  consisted  of  5977  nodes  and  11,712 
elements,  while  the  coarse  mesh  used  for  the  inverse 
calculation  had  1525  nodes  and  2928  elements.  The 
reconstructed  images  were  the  results  of  three  and 
15  iterations  for  the  method  without  and  with  the 
TVM,  respectively.  For  these  experimental  cases, 
g)<d  =  1.0  and  5  =  0.001  for  the  cases  1  and  3  and 
a>3,  =  2.0  and  <5  =  0.001  for  case  2  were  used. 

Here  we  note  that,  in  the  experimental  situation, 
the  effect  of  ultrasonic  transducer  response  should  be 
considered  in  the  PAT  image  reconstruction  because 
the  transducer  mechanical-electrical  impulse  re¬ 
sponse  as  well  as  the  transducer  aperture  effect 
may  introduce  significant  errors  in  the  forward  mod¬ 
el.  To  reduce  these  effects  in  our  calculations,  we  ap¬ 
plied  the  normalized  acoustic  pressure  distribution 
in  the  reconstruction  and  used  the  appropriate 
parameters  (initial  values,  etc.)  obtained  from  an 
optimization/searching  method  [27]. 

Figure  7  shows  the  reconstructed  absorbed  energy 
density  images  for  all  the  three  experimental  cases, 
while  Fig.  8  presents  quantitative  absorption  coeffi¬ 
cient  profiles  along  transects  through  one  target  for 
the  images  shown  in  Fig.  7.  We  see  that  considerably 
enhanced  images  are  achieved  with  the  TVM,  espe¬ 
cially  when  the  target  is  small  (case  2),  or  when  the 
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Fig.  10.  (Color  online)  Recovered  absorbed  energy  density  images  from  the  rat  brain  (a)  without  the  TVM,  (b)  with  the  TVM.  The  axes  (left 
and  bottom)  illustrate  the  spatial  scale  in  millimeters,  whereas  the  color  scale  (right)  records  the  absorbed  energy  density  in  millijoules  per 
cubed  millimeter. 


contrast  level  between  the  target  and  the  back¬ 
ground  is  low  (case  3). 

Because  there  is  no  true  or  exact  image  available 
for  the  experimental  cases,  we  used  the  contrast-to- 
noise  ratio  (CNR)  to  quantitatively  evaluate  the  ex¬ 
perimental  results.  To  compute  CNR,  we  selected 
two  regions  of  interest  (ROIs),  each  having  the  same 
size.  One  region  is  in  the  target  (referred  as  t-ROI) 
and  the  other  is  in  the  background  (denoted  as  6- 
ROI).  The  CNR  is  defined  as 


CNR  = 


I  fit)  _  fib)  I 
~a^) 


(15) 


where  and  flb] ,  respectively,  denote  the  mean  over 
the  t-ROI  and  6-ROI,  and  A7'  denotes  the  variances 
over  the  6-ROI.  We  used  solid  and  dashed  circles  to 
express  the  t-ROI  and  6-ROI  for  the  three  experimen¬ 
tal  cases.  We  calculated  the  CNR  for  each  pair  of  the 
t-  and  6-ROIs  using  Eq.  (15)  and  present  the  results 
in  Fig.  9(d).  We  see  that  the  CNR  obtained  using  the 
method  with  TVM  is  clearly  larger  than  that  using 
the  method  without  TVM,  indicating  that  the  TVM- 
based  method  is  less  sensitive  to  the  noise  effects. 

As  a  final  example,  our  methods  were  tested  using 
in  vivo  data  collected  previously  from  an  animal  (rat) 
model  of  epilepsy  [28].  Focal  seizures  were  induced 
by  microinjection  of  bicuculline  methiodide  (BMI) 
into  the  parietal  neocortex.  The  in  vivo  data  were  col¬ 
lected  using  the  same  pulsed  Nd:YAG  laser-based 
scanning  PAT  system.  In  the  reconstruction,  the  fine 
mesh  used  for  the  forward  calculation  consisted  of 
17,713  nodes  and  34,944  elements,  while  the  coarse 
mesh  used  for  the  inverse  calculation  had  4489  nodes 
and  8736  elements.  The  in  vivo  images  obtained  were 
the  results  of  two  and  five  iterations  for  the  methods 
without  and  with  the  TVM,  respectively.  We  found 
that  =  0.1  and  6  =  0.001  appeared  to  provide 
optimal  results. 

Figure  10  gives  the  recovered  images  for  the  rat 
brain  scanned  10  min  after  the  injection  of  BMI  with¬ 
out  and  with  the  use  of  the  TVM.  The  in  vivo  results 
shown  here  confirm  that  both  of  the  reconstruction 


methods  can  provide  quality  images  (the  arrow  indi¬ 
cates  the  detected  seizure  focus),  while  notable  en¬ 
hancement  in  image  quality  can  be  observed  with 
the  TVM-based  method.  For  example,  the  actual  con¬ 
tinuous  middle  cerebral  artery  (indicated  by  the 
dashed  circle)  is  shown  discontinued  in  Fig.  10(a) 
(without  the  TVM),  while  this  feature  is  clearly 
depicted  in  Fig.  10(b)  (with  the  TVM). 

4.  Conclusions 

We  have  presented  a  time-domain  FEM-based  PA 
image  reconstruction  method  that  incorporates  the 
TVM.  The  results  shown  in  this  work  indicate  that 
the  TVM-based  method  is  able  to  offer  clear  enhance¬ 
ment  in  image  reconstruction  over  the  method  with¬ 
out  the  TVM,  not  only  in  terms  of  the  location,  size, 
and  shape  of  the  target,  but  also  in  terms  of  the  ab¬ 
sorbed  energy  density  property/optical  absorption 
coefficient  values  themselves.  In  addition,  the  results 
have  shown  that  the  inclusion  of  the  TVM  in  our  re¬ 
construction  algorithm  is  highly  effective  in  the  pre¬ 
sence  of  noisy  data.  Finally,  it  is  important  to  note 
that  the  TVM-based  method  may  be  ideal  for  image 
reconstruction  involving  a  low  contrast  level  between 
the  target  and  the  background  or  small  targets  em¬ 
bedded  in  the  background. 

This  research  was  supported  in  part  by  the  Depart¬ 
ment  of  Defense  Congressionally  Directed  Medical 
Program. 
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Abstract:  Photoacoustic  tomography  (PAT)  is  an  emerging  non-invasive 
imaging  technique  with  great  potential  for  a  wide  range  of  biomedical 
imaging  applications.  However,  the  conventional  PAT  reconstruction 
algorithms  often  provide  distorted  images  with  strong  artifacts  in  cases 
when  the  signals  are  collected  from  few  measurements  or  over  an  aperture 
that  does  not  enclose  the  object.  In  this  work,  we  present  a  total-variation- 
minimization  (TVM)  enhanced  iterative  reconstruction  algorithm  that  can 
provide  excellent  photoacoustic  image  reconstruction  from  few-detector  and 
limited-angle  data.  The  enhancement  is  confirmed  and  evaluated  using 
several  phantom  experiments. 
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1.  Introduction 

Photoacoustic  tomography  (PAT)  is  an  emerging  non-invasive  imaging  technique  that 
combines  the  merits  of  high  optical  contrast  and  high  ultrasound  resolution  in  a  single 
modality  [1-3].  In  PAT,  the  Helmholtz-like  photoacoustic  wave  equation  has  been  commonly 
used  as  an  accurate  model  for  describing  laser-induced  acoustic  wave  propagation  in  tissue. 
While  a  variety  of  analytic  PAT  reconstruction  algorithms  have  been  developed  [4-6],  the 
finite  element  method  (FEM)  based  approach  appears  to  be  particularly  powerful  in  this 
regard,  because  it  can  provide  quantitative  imaging  capability  by  recovering  optical  absorption 
coefficient  [7,8],  eliminate  the  assumption  of  homogeneous  acoustic  medium  needed  in 
analytical  methods,  accommodate  object  boundary  irregularity  and  allow  appropriate 
boundary  conditions  implementations. 

A  practical  need  exists  for  reconstruction  of  photoacoustic  images  from  few 
measurements,  as  this  can  greatly  reduce  the  required  scanning  time  and  the  number  of 
ultrasound  sensors  placed  near  or  on  the  boundary  of  an  object  to  receive  the  laser-induced 
acoustic  signals.  In  addition,  in  many  practical  implementations  of  PAT  the  photoacoustic 
signals  are  recorded  over  an  aperture  that  does  not  enclose  the  object,  which  results  in  a 
limited-angle  tomographic  reconstruction  problem.  In  such  cases,  the  existing  reconstruction 
algorithms,  which  are  based  on  the  least-squares  criteria  (i.e.,  the  regularized  Newton  method) 
[9,10],  often  generate  distorted  images  with  severe  artifacts.  Total-variation-minimization 
(TVM),  on  the  other  hand,  has  proved  to  be  a  powerful  tool  for  limited-data  image 
reconstruction  in  diffraction  tomography  and  computed  tomography  (CT)  [11,12].  TVM  based 
PAT  algorithms  have  also  been  implemented  and  tested  using  numerically  simulated  limited- 
view  data  [13-16].  Here  we  describe  an  unconstrained  TVM  FEM-based  iterative 
reconstruction  algorithm,  and  for  the  first  time  present  experimental  evidence  that  the  TVM- 
based  algorithm  offers  excellent  photoacoustic  image  reconstruction  from  few-detector  and 
limited-angle  data. 

This  paper  is  organized  as  follows.  In  Section  2,  the  FEM  based  PAT  reconstruction 
algorithm  and  the  unconstrained  TVM  scheme  are  reviewed  briefly.  The  experimental 
validation  of  our  TVM  FEM-based  algorithm  is  presented  in  Section  3.  The  conclusions  are 
made  in  Section  4. 

2.  Method 

We  first  briefly  describe  the  existing  FEM-based  photoacoustic  reconstruction  algorithm 
detailed  elsewhere  [17].  The  time-domain  photoacoustic  wave  equation  in  tissue  can  be 
described  as  follows: 


v  ’  Vq  dt2  Cpdt 

where  p  is  the  pressure  wave;  v0  is  the  speed  of  acoustic  wave  in  the  medium;  //  is  the  thermal 
expansion  coefficient;  Cp  is  the  specific  heat;  <D  is  the  absorbed  energy  density; 
J(t)  =  S(t-t0)  is  assumed  in  our  study. 
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To  form  an  image  from  a  presumably  uniform  initial  guess  of  the  absorbed  energy  density 
distribution,  we  need  a  method  of  updating  ®,  the  desired  quantity  from  its  starting  value.  This 
update  is  accomplished  through  the  least-squares  minimization  of  the  following  functional: 

F(p’®)=i{p0j-p:) >  p) 

j= i 

where  p°  and  pCj  are  the  measured  and  computed  acoustic  field  data  for  i  =  \,2,---M 

boundary  locations.  Using  the  regularized  Newton  method,  we  obtained  the  following  matrix 
equation  for  updating  O: 

(3r3  +  Al)AZ  =  3r  (//-//),  (3) 

where  p°  =  [pi,P2,'"Pm)  and  P°  =  Pm  )  >  A/  >s  the  update  vector  for  the 

absorbed  optical  energy  density;  3  is  the  Jacobian  matrix  formed  by  dpjdt D  at  the  boundary 
measurement  sites;  1  is  the  regularization  parameter  determined  by  combined  Marquardt  and 
Tikhonov  regularization  schemes;  and  I  is  the  identity  matrix. 

We  now  incorporate  the  total  variation  of  <5  as  a  penalty  term  by  defining  a  new  functional 
[18]: 

#(/>,©)  =  F(/»,©)  +  £(©).  (4) 

Here  Z,(®)  =  |V<t>|‘  +  S2  dxdy  is  the  penalty  term,  and  <»„,  and  S  are  typically  positive 

parameters  that  need  to  be  determined  numerically.  The  minimization  of  Eq.  (4)  can  be 
realized  by  the  differentiation  of  F  with  respect  to  each  nodal  parameter  that  constitutes  the  ® 
distribution  and  by  setting  each  of  the  resulting  expression  to  zero,  leading  to  the  following 
system  of  equations: 


8F 

5®,. 


M  r)T)C 

o  (<=u-n 


(5) 


where  Vj  =  0Z/5®,.  .  Again,  through  the  regularized  Newton  method,  the  following  matrix 
equation  for  TVM-based  inversion  is  obtained  [18]: 

(5t5  +  R  +  M)aX  =  3t(p°-p')-V,  (6) 


where  V  is  formed  by  cL / f®  and  R  is  formed  by  dV / 1®  .  At  this  point,  our  solution 
procedure  for  solving  Eq.  (6)  and  the  regularization  parameter  selection  are  identical  to  those 
described  previously.  Hence,  it  becomes  clear  that  the  only  additions  to  our  new  algorithm 
result  from  the  assembly  of  matrix  R  and  the  construction  of  column  vector  V. 

3.  Results 

In  this  section  our  TVM  enhanced  reconstruction  algorithm  is  tested  and  evaluated  using 
phantom  experimental  data.  For  comparative  purposes,  reconstruction  results  without  the 
TVM  enhancement  are  also  presented. 

The  experimental  setup  used  for  collecting  the  phantom  data  was  a  pulsed  ND:  YAG  laser 
based  single  transducer  (1MHz)  scanning  system,  which  was  described  in  detail  elsewhere  [8], 
Three  phantom  experiments  were  conducted.  In  the  first  two  experiments,  we  embedded  one 
or  two  objects  with  a  size  ranging  from  3  to  0.5  mm  in  a  50  mm-diameter  solid  cylindrical 
phantom.  The  absorption  coefficient  of  the  background  phantom  was  0.01  mm-1,  while  the 
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absorption  coefficient  of  the  target(s)  was  0.03  mm-1.  In  the  last  experiment,  we  used  a  single¬ 
target-containing  phantom,  aiming  to  test  the  capability  of  detecting  target  having  low  optical 
contrasts  relative  to  the  background  phantom.  In  this  case,  the  target  had  an  absorption 
coefficient  of  0.015  mm"1.  The  reduced  scattering  coefficients  of  the  background  phantom  and 
targets  were  1.0  and  3.0  mm-1  for  the  first  two  experiments,  and  1.0  and  2.0  mm-1  for  the  last 
experiment.  In  the  reconstructions,  we  used  a  dual-meshing  scheme  [19]  where  the  fine  mesh 
used  for  the  forward  calculation  consisted  of  5977  nodes  and  11712  elements,  while  the 
coarse  mesh  used  for  the  inverse  calculation  had  1525  nodes  and  2928  elements.  All  the 
images  obtained  from  the  method  without  the  TVM  are  the  results  of  three  iterations,  while 
those  obtained  from  the  TVM-enhanced  method  are  the  results  of  twenty  or  more  iterations. 
Parallel  code  was  used  to  speed  up  these  calculations  on  Beowulf  clusters  with  8  CPUs.  The 
computational  speed  can  be  further  increased  if  clusters  with  more  than  8  CPUs  are  used.  As 
mentioned  above,  the  parameters  a>0  and  8  were  determined  through  numerical 
experimentation.  From  our  experience,  a  constant  value  of  8  =  0.001  was  sufficient  for  the 
current  experimental  studies,  while  the  value  of  co,h  is  related  to  the  signal-to-noise  ratio  of  the 
measurements.  For  our  experimental  cases,  co^  =1.0  and  8  =  0.001  for  the  first  case,  and 
cOq  =  0.5  and  8  =  0.001  for  the  second  and  third  cases  were  used. 

The  reconstruction  results  based  on  few-detector  data  from  the  three  experimental  cases 
are  shown  in  Figs.  1,  2  and  3,  respectively.  In  each  figure,  the  images  in  the  top  and  bottom 
rows,  respectively,  present  the  recovered  absorbed  energy  density  images  using  the  algorithm 
without  and  with  the  TVM  where  120,  60,  30,  and  15  detectors  were  equally  distributed  along 
the  surface  of  the  circular  background  region.  As  expected,  the  conventional  algorithm 
without  the  TVM  provides  high  quality  images  only  when  the  number  of  detectors  was 
relatively  sufficient  (Figs.  la,b,  Figs.  2a, b,  and  Figs.  3a, b).  The  images  reconstructed  from 
few-detector  data  by  the  conventional  algorithm  without  the  TVM  contained  severe  artifacts 
and  distortions  (Figs.  lc,d.  Figs.  2c, d,  and  Figs.  3c,d).  Considerably  enhanced  images  are 
achieved  using  the  method  with  the  TVM,  especially  when  the  number  of  the  detectors  is 
reduced  to  30  or  15  (Figs.  lg,h.  Figs.  2g,h,  and  Figs.  3g,h).  We  also  note  that  the 
computational  efficiency  of  our  TVM  based  algorithm  for  recovering  a  small  absorber  (Fig.  2) 
and  a  larger  absorber  (Fig.  3)  is  similar  when  the  number  of  the  detectors  is  insufficient. 

Images  reconstructed  from  the  limited-angle  data  for  the  first  case  using  the  two  methods 
are  displayed  in  Fig.  4  where  the  top  and  bottom  rows,  respectively,  show  the  recovered 
absorbed  energy  density  images  using  the  algorithm  without  and  with  the  TVM  when  120 
detectors  over  360°,  60  detectors  over  180°,  and  30  detectors  over  90°,  were  equally 


Fig.  1.  Reconstructed  photoacoustic  images  based  on  few-detector  data  for  case  1.  (a),  120 
detectors,  without  TVM.  (b),  60  detectors,  without  TVM.  (c),  30  detectors,  without  TVM.  (d), 
15  detectors,  without  TVM.  (e),  120  detectors,  with  TVM.  (f),  60  detectors,  with  TVM.  (g),  30 
detectors,  with  TVM.  (h),  15  detectors,  with  TVM. 
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Fig.  2.  Reconstructed  photoacoustic  images  based  on  few-detector  data  for  case  2.  (a),  120 
detectors,  without  TVM.  (b),  60  detectors,  without  TVM.  (c),  30  detectors,  without  TVM.  (d), 
15  detectors,  without  TVM.  (e),  120  detectors,  with  TVM.  (f),  60  detectors,  with  TVM.  (g),  30 
detectors,  with  TVM.  (h),  15  detectors,  with  TVM. 


Fig.  3.  Reconstructed  photoacoustic  images  based  on  few-detector  data  for  case  3.  (a),  120 
detectors,  without  TVM.  (b),  60  detectors,  without  TVM.  (c),  30  detectors,  without  TVM.  (d), 
15  detectors,  without  TVM.  (e),  120  detectors,  with  TVM.  (f),  60  detectors,  with  TVM.  (g),  30 
detectors,  with  TVM.  (h),  15  detectors,  with  TVM. 


Fig.  4.  Reconstructed  photoacoustic  images  based  on  limited-angle  data  for  case  1.  (a),  120 
detectors  over  360°,  without  TVM.  (b),  60  detectors  over  180°,  without  TVM.  (c),  30  detectors 
over  90°,  without  TVM.  (d),  120  detectors  over  360°,  with  TVM.  (e),  60  detectors  over  180°, 
with  TVM.  (f),  30  detectors  over  90°,  with  TVM. 

distributed  along  the  surface  of  the  circular  background  region.  Again,  strong  artifacts  and 
distortions  exist  in  the  images  recovered  with  the  method  without  the  TVM  when  the  data 
collected  is  angle-limited  (Figs.  4b, c),  while  considerably  improved  images  are  reconstructed 
using  the  TVM  enhanced  algorithm  (Figs.  4e,f). 
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In  order  to  evaluate  quantitatively  the  reconstruction  quality  using  the  method  with  and 
without  TVM  enhancement,  we  use  the  following  universal  quality  index  (UQI)  [20]: 


UQl{  f',f0} 


2Cov{f‘,f0}  2  f'f° 


(7) 


Here  the  image  f  can  be  interpreted  as  vectors  of  size  N:  f  =  (f1,f2,'"fN)T ,  where  N  denotes 
the  number  of  image  data  acquired  from  the  FEM  based  algorithm.  fJ  is  the  image  means, 
aJ  is  the  variances,  and  Covjf'jf0]  is  the  covariances  over  the  whole  image  domain,  where y 

=  0  and  1.  UQI  measures  the  image  similarity  between  the  reconstructed  (f1)  and  reference 
( f° )  images,  and  its  value  ranges  between  0  and  1.  The  value  of  the  UQI  is  closer  to  1  when 
the  reconstructed  image  is  more  similar  to  the  reference  image. 

We  calculated  the  UQIs  for  the  three  cases  presented  above,  and  the  results  are  given  in 
Fig.  5  where  Figs.  5a  and  5b  show  the  results  from  the  few-detector  data  based  on  the  method 
without  and  with  the  TVM,  respectively,  while  Fig.  5c  presents  the  results  from  the  limited- 
angle  data  for  the  first  case  based  on  the  two  methods.  Here  the  image  data  recovered  by  120 
detectors  was  regarded  as  the  reference  one.  The  computed  UQIs  shown  in  Fig.  5  confirm  the 
observation  that  the  TVM  enhanced  PAT  algorithm  provides  significantly  better  image  quality 
compared  to  the  conventional  method  without  the  TVM. 


(a)  0>)  (c) 

Fig.  5.  UQIs  calculated  from  the  recovered  images  for  the  three  cases  based  on  few-detector 
data  without  TVM  (a)  and  with  TVM  (b),  and  for  case  1  based  on  limited  angle  data  with  and 
without  TVM. 


4.  Conclusions 


In  this  work,  we  have  implemented  and  evaluated  an  unconstrained,  total-variation- 
minimization  method  for  time-domain  FEM  based  PAT.  This  study  has  confirmed  that  in  the 
situation  that  the  data  is  collected  from  few  measurements  or  over  an  aperture  that  does  not 
enclose  the  object,  the  developed  TVM  based  PAT  algorithm  provides  considerably  improved 
image  reconstruction  compared  to  the  conventional  methods.  The  application  of  this  TVM 
enhanced  reconstruction  algorithm  will  significantly  reduce  the  number  of  ultrasound 
transducers  and  scanning  time  needed  for  high  quality  photoacoustic  image  reconstruction. 
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ABSTRACT 


While  brain  imaging  and  electrophysiology  independently  play  a  central  role  in  today’s 
neuroscience  and  in  the  evaluation  of  neurological  disorders,  a  single  noninvasive  modality 
that  offers  both  high  spatial  and  temporal  resolution  at  centimeter  scale  depths  is  currently 
not  available.  Here  we  show  in  an  acute  epilepsy  rat  model  that  photoacoustic  tomography 
(PAT)  can  noninvasively  track  seizure  brain  dynamics  with  both  high  spatial  and  temporal 
resolution,  and  at  a  depth  that  is  clinically  relevant.  The  noninvasive  yet  whole  surface  and 
depth  capabilities  of  the  photoacoustic  imaging  system  allowed  for  the  first  time  to  actually 
see  what  is  happening  during  ictogenesis  in  terms  of  seizure  onset  and  spread.  Both  seizure 
onset  and  propagation  were  tomographically  detected  at  a  spatial  resolution  of  150  pm  and  a 
temporal  resolution  of  300  ms,  respectively.  The  networks  associated  with  the  epileptiform 
events  were  identified  using  Granger  Causality,  which  lend  support  to  the  theory  that 
seizure  onset  and  spread  involves  a  rich  interplay  between  multiple  cortical  and  subcortical 
foci  during  the  onset  and  spread  of  focal  epilepsy.  Dynamical  changes  of  vasculature  during 
epileptiform  events  were  also  detected  with  high  spatiotemporal  resolution.  Together,  these 
findings  suggest  that  PAT  represents  a  powerful  tool  for  noninvasively  mapping  seizure  onset 
and  propagation  patterns,  and  the  ‘functional’  connectivity  within  epileptic  brain  networks. 

Keywords:  neuroimaging;  epilepsy;  localization;  neural  networks;  microvasculature 
Abbreviations:  PAT=photoacoustic  tomography;  EEG=  electroencephalography; 

BMI=  Bicuculline  methiodide  influences;  ROI=region  of  interest 
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INTRODUCTION 


Epilepsy  is  a  common,  chronic,  neurological  disorder  characterized  by  seizures.  Three 
percent  of  people  will  be  diagnosed  with  epilepsy  at  some  time  in  their  lives  (Hauser  and 
Kurland  1975,  Hauser  1992).  Indeed,  approximately  50  million  people  worldwide  have 
epilepsy,  and  20  to  30  percent  of  these  patients  are  refractory  to  all  forms  of  medical 
treatment  (Hauser  1 992)  .  Seizures  are  transient  signs  and  symptoms  of  abnormal,  excessive, 
or  hypersynchronous  neuronal  activity  in  the  brain  (Fisher  et  al.  2005).  In  most  cases, 
seizures  are  controlled,  although  not  cured,  with  anticonvulsant  medication.  Seizure  types 
are  classified  according  to  whether  the  source  of  the  seizure  is  localized  (partial  or  focal  onset 
seizures)  or  distributed  (generalized  seizures)  (Greenfield  2011).  Localization-related 
epilepsies  arise  from  an  epileptic  focus.  A  partial  seizure  may  spread  within  the  brain,  a 
process  known  as  secondary  generalization.  Generalized  epilepsies,  in  contrast,  arise  from 
many  independent  foci  (multifocal  epilepsies)  or  from  epileptic  circuits  that  involve  the 
whole  brain.  In  epilepsies  of  unknown  localization,  it  remains  unclear  as  to  whether  they 
arise  from  a  portion  of  the  brain  or  from  more  widespread  circuits. 

For  those  patients  with  medically  intractable  focal  epilepsy,  the  best  treatment  option  is 
resective  brain  surgery  (Birbeck  et  al.  2002,  Duncan  et  al.  2006,  Berg  et  al.  2007).  Although 
several  factors  can  impact  the  success  of  epilepsy  surgery,  the  primary  reason  of  failure  is  the 
incomplete  mapping  of  the  local  epilepsy  network  which  results  in  incomplete  resection  of 
epileptogenic  foci  (Engel  2004,  Jeha  et  al.  2007).  Much  of  the  attention  on  epilepsy  surgery 
has  been  directed  at  identifying  single  neuronal  populations.  This  approach  has,  in  many 
cases,  led  to  failed  surgical  outcomes,  because  seizures  typically  involve  groups  of  neurons 
interacting  both  locally  and  across  several  cortical  and  subcortical  brain  regions.  A  better 
understanding  of  the  regional  interactions  occurring  at  the  site  of  seizure  onset  and  spread 
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may  provide  important  insights  about  the  pathophysiology  of  seizures  and  aid  with  accurate 
brain  mapping  and  resection  of  the  epileptic  focus. 

The  epileptic  focus  is  defined  as  the  portion  of  the  brain  that  serves  as  the  irritant, 
driving  the  epileptic  response.  The  seizure  focus  is  generally  considered  as  the  place  where  a 
seizure  starts  and  the  target  for  surgical  intervention.  The  common  belief  is  that  the  seizure 
starts  the  same  way,  and  that  the  entire  focus  acts  in  unison.  In  theory,  removing  the  focus 
should  result  in  a  patient’s  seizures  being  cured.  However,  there  is  much  evidence  to  suggest 
that  the  focus  is  more  of  a  region  of  seizure  onset  with  a  number  of  sites  that  can  act 
independently  to  initiate  seizures  (Thom  et  al.  2010).  Seizures  in  animal  models  and  in 
people  often  have  a  multifocal  or  broadly  synchronized  onset.  The  best  evidence  for 
multifocality  within  the  seizure  onset  zone  comes  from  surgical  experience  with  intracranial 
monitoring.  The  challenge  of  mapping  the  epileptic  focus  stems  from  the  observation  that  the 
pathology  associated  with  focal  epilepsy  is  often  distributed  across  a  number  of  brain  sites 
(Bertram  2009)  and  that  current  diagnostics  methods  frequently  fall  short  of  identifying  such 
sites.  Animal  studies  indicate  that  the  neurons  involved  in  the  epileptic  circuitry  have 
enhanced  excitability  throughout  (Bertram  et  al.  1998,  Fountain  et  al.  1998,  Mangan  et  al. 
2000).  The  implication  of  these  observations  is  that  each  of  the  sites  could  act  independently 
to  initiate  a  seizure  or,  potentially,  to  drive  another  site  into  a  seizure.  Thus,  in  focal  epilepsy, 
one  may  view  a  cortical  region  as  a  broad  seizure  onset  zone,  with  the  potential  that  multiple 
foci  can  act  as  a  seizure  focus  for  any  given  seizure. 

Much  of  our  understanding  about  focal  seizure  circuitry  comes  from 
electrophysiological  recording  methods.  Although  electrophysiology  is  currently  the  ‘gold 
standard’  in  mapping  the  epileptic  focus,  it  is  often  inadequate  to  define  the  boundary  of  the 
epilepsy  circuitry  due  to  spatial  sampling  limitations  and  volume  conduction.  What  we  truly 

desire  is  a  brain  mapping  modality,  which  could  give  a  high-resolution  real  time  spatial  and 
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temporal  ‘read  out’  of  the  dynamics  of  cortical  processing  and  seizures.  It  has  been  well 
established  that  optical  contrast  is  highly  sensitive  to  neuronal  activity  (Hill  and  Keynes 
1949,  Grinvald  et  al.  1988).  The  high  optical  contrast  is  largely  due  to  the  changes  in  blood 
volume  and  blood  oxygenation.  When  used  intraoperatively,  optical  imaging  or  intrinsic 
optical  signal  provided  excellent  surface  maps  of  epileptic  focus  (Haglund  et  al.  1992, 
Haglund  and  Hochman  2004).  Near-infrared  spectroscopy  has  also  been  used  to  obtain 
cerebral  blood  volume  and  blood  oxygenation  in  epilepsy  patients  (Steinhoff  et  al.  1996, 
Sokol  et  al.  2000,  Watanabe  et  al.  2000).  However,  the  major  limitation  of  Near-infrared 
spectroscopy  and  or  intrinsic  optical  signa  is  that  each  provides  only  surface  depth 
information.  Greater  depth  information  can  be  obtained  by  using  tomographic  reconstruction 
methods  such  as  diffuse  optical  tomography  (Oleary  et  al.  1995,  Jiang  et  al.  1996,  Bluestone 
et  al.  2001,  Boas  et  al.  2001).  Whereas  diffuse  optical  tomography  has  low  spatial  resolution, 
laser-induced  photoacoustic  tomography  (PAT)  has  both  superior  spatial  and  temporal 
resolution.  PAT  detects  absorbed  photons  ultrasonically  by  employing  the  photoacoustic 
effect.  It  combines  both  high  contrast  and  spectroscopic  specificity  based  on  the  optical 
absorption  of  both  oxy-  and  deoxy-hemoglobin  with  high  ultrasonic  spatial  resolution 
(Kruger  and  Liu  1994,  Laufer  et  al.  2007,  Xiang  et  al.  2007,  Yuan  et  al.  2007).  Relative  to 
other  optical  imaging  modalities,  PAT  has  the  advantage  of  mitigating  both  scalp  and  skull 
light  scattering  by  a  factor  of  -1,000.  The  end  result  is  that  PAT  allows  for  high  spatial 
resolution  imaging  of  brain  at  a  depth  considerably  beyond  the  soft  depth  limit  of 
conventional  optical  imaging  techniques  such  as  confocal  microscopy  (Sipkins  et  al.  2005), 
two-photon  microscopy  (Denk  et  al.  1990),  and  optical  coherence  tomography  (Huang  et  al. 
1991). 

In  this  study,  PAT  was  employed  to  image  seizures  in  an  experimental  acute  bicuculline 

methiodide  model  of  focal  epilepsy.  Bicuculline  is  a  light-sensitive  competitive  antagonist  of 
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GABAa  receptors  that  mimics  focal  epilepsy  when  applied  to  brain  tissue.  During  focal 
application  of  bicuculline  into  the  brain  cortex,  brains  were  imaged  noninvasively  with  a 
novel  PAT  system  that  has  three  orders  of  magnitude  higher  temporal  resolution  and  four¬ 
fold  higher  spatial  resolution  relative  to  our  previous  PAT  prototype(Zhang  et  al.  2008).  The 
high  spatiotemporal  properties  of  this  PAT  imaging  tool  allowed  for  identification  of  both 
bicuculline  elicited  focal  seizure  onset  and  spread.  Off-line,  we  employed  measures  of  brain 
connectivity  to  further  identify  the  functional  anatomy  implicated  in  focal  cortical  seizures. 
The  high  spatial  and  temporal  sampling  of  the  novel  PAT  system  allowed  for  the  first  time 
the  complete  mapping  of  an  epileptiform  event  in  vivo.  It  is  also  the  first  report  of  mapping 
an  epileptiform  event  at  depths  well  below  the  cortex.  These  findings  build  onto  the 
collection  of  advantages  of  noninvasive  PAT  in  brain  dynamics  and  epilepsy.  Information 
from  this  study  can  be  applied  to  improve  our  understanding  of  epileptic  networks  and 
mapping  in  epilepsy.  In  terms  of  impact,  this  is  the  first  demonstration  of  a  potentially 
clinical  useful  adjunct  in  the  clinical  surgical  evaluation  of  cortical  epilepsy. 

MATERIALS  AND  METHODS 

Animals 

Male  Sprague-Dawley  rats  (Harlan  Labs,  Indianapolis,  IN)  weighing  50-60  g  on  arrival 
were  allowed  one  week  to  acclimate  to  the  12-h  light/dark  cycle  and  given  food  and  water  ad 
libitum.  All  procedures  were  approved  by  the  University  of  Florida  Animal  Care  and  Use 
Committee  and  conducted  in  accordance  with  the  National  Institutes  of  Health  Guide  for  the 
Care  and  Use  of  Experimental  Animals. 

Electrode  implantation  surgery 

Animals  were  anesthetized  intraperitoneally  with  1  g/kg  of  body  weight  dose  of 

urethane.  Two  300  pm  diameter  stainless  steel  screw  electrodes  were  implanted  within  the 
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skull  for  obtaining  subdural  multichannel  cortical  local  field  potentials  data  (-0.3  mm 
posterior,  3  mm  lateral  (right)  of  bregma,  1  mm  ventral)  based  on  coordinates  from  a  rat  brain 
atlas(Paxinos  1998).  One,  300  pm  diameter  stainless  steel  screw  electrode  (FHC,  Bowdoin, 
ME)  was  implanted  as  a  reference  electrode  into  the  midline  occipital  bone.  Cortical  local 
field  potentials  were  obtained  using  a  Tucker  Davis  Pentusa  (Tucker  Davis  Technologies, 
Alachua,  FL)  neural  recording  system  at  12  kHz,  digitized  with  16  bits  of  resolution,  and 
band  pass  filtered  from  0.5  to  6  kHz. 

Induction  of  seizures 

Rats  (n=10)  received  10  pi  of  1.9  mM  bicuculline  methiodide  influences  (BMI)  and 
10  pi  normal  saline  into  the  left  and  right  parietal  cortex,  respectively.  The  infusion  was 
performed  through  the  previously  implanted  electrode  sites  at  a  rate  of  0.3  pL/min.  The 
infusion  system  consisted  of  a  100  pL  gas-tight  syringe  (Hamilton,  Reno,  NV)  driven  by  a 
syringe  pump  (Cole-Parmer,  Vernon  Hills,  IL)  connected  to  polyaryletheretherketone 
(PEEK)  tubing  (ID  =  0.381mm,  OD=  0.794mm,  length<0.5  m,  Upchurch  Scientific,  Oak 
Harbor,  WA).  The  PEEK  tubing  was  coupled  to  a  silica  cannula  (ID  =  50  pm,  OD=  147  pm, 
Polymicro  Technologies,  Phoenix,  AZ)  via  a  microfluidic  connector.  Cortical  local  field 
potentials  were  recorded  5  min  before  each  injection  and  continued  for  up  to  30  min 
thereafter. 

Image  analysis 

A  custom  software  utility  was  written  and  incorporated  into  the  computer  software 
package  MATLAB  (MathWorks,  Massachusetts,  USA)  to  analyze  and  display  recorded  data. 
This  software  enabled  the  reconstruction  of  the  PAT  images  and  the  determination  of  the 
blood  vessel  diameter.  Amira  (version  5.3.3,  TGS  Template  Graphics  Software)  was  used  for 
three-dimensional  reconstruction  of  the  seizure  foci.  In  connectivity  analysis,  the  time  series 
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data  were  selected  from  200  sets  of  PAT  images  sampled  at  3Hz  within  the  seizure  onset 
period,  and  at  each  time  point,  we  chose  9  (regions  of  interest)  ROI  ( R I  ~R9)  and  then  the  PA 
signal  was  averaged  within  each  ROI,  as  shown  in  Fig.  4ci.  Changes  in  photoacoustic  signals 
were  quantified  as  -AA/A  in  Fig.  4b,  right.  The  PA  data  shown  in  Figs.  2-5  were  not 
averaged.  The  image  color  scale  was  determined  by  the  photoacoustic  signal  intensity  in 
arbitrary  units  (between  0  to  1). 

RESULTS 

PAT  imaging 

Light  from  a  Ti:  Sapphire  laser  tunable  (690  to  950  nm)  was  delivered  through  the  skull 
to  the  brain  through  an  optical  fiber  (Fig.  la).  The  energy  of  each  laser  pulse  was  detected  by 
a  photodiode  for  calibration.  A  1 92  element  full-ring  transducer  array  was  used  to  capture  the 
photoacoustic  (PA)  signals  generated  by  the  laser  light.  The  192  channel  data  acquisition 
system  consisted  of  preamplifiers,  secondary  stage  amplifiers  (for  optimizing  the  signal-to- 
noise  ratio),  and  a  3:1  electronic  multiplexer  coupled  with  a  64-channel  analog-to-digital 
converter.  The  estimated  PAT  data  acquisition  system  speed  was  0.33  s/frame.  Data 
acquisition  speed  was  limited  in  part  by  the  10  Hz  laser  repetition  (supplemental  figure  1, 
Movie  S3).  The  spatial  resolution  of  the  PAT  system  was  inversely  related  to  the  bandwidth 
of  the  ultrasonic  transducer.  Each  ultrasonic  detector  had  a  5-MHz  central  frequency  and  a 
70%  nominal  bandwidth  with  a  diameter  of  6mm  (Blatek,  Inc.,  PA,  U.S.A.).  The 
noninvasive  PAT  image  of  the  rat  cortical  vasculature  (Fig.  lb)  matched  well  with  an 
anatomical  photograph  (Fig.  lc)  obtained  after  the  PA  imaging.  The  brain  structures 
including  the  middle  cerebral  artery,  right  hemispheres,  left  hemispheres,  left  olfactory 
bulbs,  and  right  olfactory  bulbs  are  clearly  shown  in  the  PAT  image.  Numbers  1-5 
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marked  in  the  image  indicate  that  the  micro-blood  vessels  with  a  diameter  of  less  than 
100pm  are  also  seen,  which  again  correspond  well  with  the  rat  brain  photograph.  This 
allowed  for  imaging  objects  of  50pm  in  diameter  with  a  spatial  resolution  of  150pm 
(supplemental  figure  2). 

Thirty  minutes  following  the  microelectrode  placement,  PAT  imaging  was  performed  to 
generate  baseline  tissue  absorption  maps,  visualize  micro  blood  vessels,  and  morphological 
cortical  landmarks  (Figs  2b,  c).  Subsequently,  rats  (n=10)  received  10  pi  of  1.9  mM 
bicuculline  methiodide  and  10  pi  normal  saline  into  the  left  [-1.0AP,  2.5ML,  2.4DV]  and 
right  [-1.0AP,  2.5  ML,  0.4DV]  parietal  cortex,  respectively,  based  on  coordinates  from  a  rat 
brain  atlas  (Paxinos  1998).  Immediately  following  the  focal  infusions,  PAT  imaging  was 
repeated  to  visualize  changes  in  the  injection  site  tissue  absorption.  An  increase  in  tissue 
absorption  was  observed  in  the  bicuculline  methiodide  cortical  injection  site  and  surrounding 
region  (Fig.  2c,  left),  but  not  in  the  saline  injection  site  or  surrounding  region  (Fig.  2c,  right). 
PAT  scanning  for  subcortical  changes  was  also  performed.  Slices  5,  7,  and  9  (Fig.  2d)  are 
the  images  obtained  3,  5  and  7  mm  below  the  scalp,  respectively.  Fig.  2e  is  the  three 
dimensional  rendering  of  the  epileptic  foci  obtained  from  different  tomographic  layers 
(movie  SI).  Each  image  exhibited  the  variable  patterns  of  seizure  onset  and  propagation  both 
at  the  cortical  and  subcortical  regions.  Seizures  were  confirmed  by  concomitant  time-locked 
PAT/video-electroencephalography  (Fig.  2a).  Experiments  were  repeated  for  each  animal  at 
2-hour  intervals. 

Real-time  monitoring  of  epileptic  events 

Epileptiform  events  were  recorded  with  PAT  from  a  focal  region  of  interest  of  -2x3 
mm  (Fig.  3a).  We  observed  a  significant  optical  absorption  change  directly  associated  with 
de-oxygenation  at  a  wavelength  of  755  nm  (Fig.  3c).  Furthennore,  at  the  seizure  onset  at 
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Imin  6.267  s,  the  seizure  focus  measured  0.2053  mnr  and  increased  to  0.5004  mnr  as  the 
electrographic  seizure  time-series  increased  in  frequency  and  amplitude  (Fig.  3b). 
Corresponding  rate  changes  (0.33  s)  spike  and  wave  discharges  and  PA  images  were 
observed  in  each  experimental  trial.  The  seizure  onset  dynamics  captured  photoacoustically 
are  best  appreciated  in  movies  displaying  how  the  seizure  was  generated  and  varied  over  time 
(movie  S2).  PAT  images  suggest  that  in  addition  to  ictal  spread  from  the  primary  focus, 
homotopic  foci  are  seen  in  the  contralateral  parietal  cortex  (Fig.  4ai_4).  The  PA  signal  in  the 
contralateral  foci  was  smaller  in  magnitude  and  delayed  in  time  compared  to  the  signal 
recorded  from  the  primary  seizure  foci.  We  also  observed  a  region  of  inverted  optical  signal 
immediately  surrounding  the  seizure  onset  zone  (Fig  4b,  left).  We  performed  analysis  of  the 
optical  signal  at  a  distance  of  at  least  2  mm  away  from  the  edge  of  the  focus.  Our  results 
demonstrate  that  the  inverse  optical  signal  was  inversely  related  to  the  optical  signal  recorded 
from  the  focus  (Fig.  4b,  right). 

To  evaluate  the  dynamic  interactions  within  the  seizure  circuitry  and  capture  the 
associated  networks  we  used  the  frequency  domain  Granger  causality  (Brovelli  et  al.  2004, 
Granger  1969,  Wang  et  al.  2007) .  PAT  time  series  data  was  selected  from  200  sets  of  PAT 
images  sampled  at  3  Hz  within  the  seizure  onset  period.  At  each  time  point  we  chose  9  ROIs 
(Fig.  4cj)  and  used  the  averaged  optical  absorption  within  each  region  for  the  network 
analysis.  Nine  sets  of  time  series  analysis  were  classified  into  3  groups  including  the  primary 
ictal  onset  area  and  corresponding  4  surrounding  regions  of  interest  including  the  primary  or 
initiating  seizure  focus,  contralateral  homotopic  foci,  and  2  regions  of  interest  surrounding 
the  primary  focus  in  the  contralateral  cortex.  The  middle  cerebral  artery,  primary  seizure 
focus,  and  secondary  homotopic  seizure  focus  were  also  analyzed.  Results  from  group  1 
analysis  demonstrate  the  causal  influence  from  the  primary  focus  to  the  surrounding  regions 


of  interest  (Fig.  4b,  left).  The  negative  zero-lag  correlation  (-0.68)  between  both  the  primary 
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and  secondary  foci  showed  that  neural  activities  within  these  two  regions  changed  in  opposite 
direction  over  time.  We  postulate  that  the  inverted  change  in  optical  signal  was  likely  due  to 
the  opposite  oxygen  delivery  or  blood  flow(Schwartz  and  Bonhoeffer  2001).  Figure  4c2 
shows  the  Granger  causality  analysis  for  group  2  data  where  we  see  that  the  primary  and 
secondary  contralateral  brain  foci  influenced  each  other  mutually  with  a  stronger  influence  of 
the  primary  focus  relative  to  the  contralateral  focus.  Moreover,  we  suggest  that  the  primary 
seizure  focus  may  influence  the  contralateral  foci  as  previously  suggested  using  optical 
intrinsic  imaging  in  epilepsy  (Khalilov  et  al.  2003)  and  microelectrode  recordings  of 
hippocampus  local  field  potentials(Zhou  et  al.  2009,  Cadotte  et  al.  2010).  Finally,  from 
group  3  analyses  (Fig.  4C3)  we  found  that  the  primary  focus  strongly  influenced  the 
hemodynamic  change  in  the  middle  cerebral  artery,  which  in  turn  showed  influence  on  the 
contralateral  homo  topic  foci. 

Dynamical  changes  of  vasculature  during  interictal  discharges 

Fig.  5a  shows  image  of  the  rat  cortex  vasculature  through  the  intact  scalp  and  skull.  The 
red  arrow  indicates  the  microvascularity  along  the  direction  marked  by  the  yellow  dashed 
line,  representing  a  cross-sectional  scanning  using  a  50  MHz  ultrasound  transducer.  Positive 
and  negative  acoustic  peaks  induced  by  a  25  pm  blood  vessel  were  sorted  out  as  target 
signals  to  track  the  change  of  the  vessel  diameter.  Typical  PA  signals  from  the  targeted  vessel 
at  different  times  show  apparent  vessel  vasomotion  by  ~40%  during  the  interictal  discharges 
(Fig.  5b).  Local  field  potential  recordings  (Fig.  5c)  showed  that  interictal  spikes  had  a  strong 
correlation  with  the  discrete  change  in  vessel  size  observed  photoacoustically  (Fig.  5d;  error 
bars  ±  s.d.  were  calculated  from  20  consecutive  spikes). 
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DISCUSSION 


The  main  finding  is  that  PAT  imaging  is  a  powerful  tool  for  noninvasively  mapping  seizure 
dynamics  with  both  high  spatial  and  temporal  resolution.  This  study  investigated  the 
hemodynamic  changes  during  focal  seizure  onset  and  propagation  in  an  acute  model  of  focal 
epilepsy.  The  experimental  results  suggest  that  the  increase  in  local  and  surrounding  brain 
tissue  absorption  was  due  to  the  focal  bicuculline  methiodide  induced  seizure  rather  than 
other  factors  such  as  stab  cortex  tissue  wound  injection.  Epileptiform  events,  including 
interictal  spikes,  ictal  onset  and  spread  were  identified  in  near  real  time.  To  the  best  of  our 
knowledge,  this  is  the  first  experimental  evidence  for  millisecond  temporal  resolution, 
micron  scale  imaging,  and  centimeters  depth,  of  focal  seizure  onset  and  spread  by  a 
noninvasive  technique. 

Quantitative  photoacoustic  image  reconstruction  can  be  used  to  resolve  the  light 
scattering  problem  with  the  deeper  penetration  in  the  tissue  based  on  compensating  the 
reconstructed  image  with  a  pre-calculated  optical  fluence  distribution(Yuan  et  al.  2007,  Yao 
et  al.  2009).  It  is  possible  to  achieve  higher  resolution  at  cellular  and  even  sub-cellular  levels 
to  image  the  neuronal  circuitry  in  vivo  by  use  of  photoacoustic  microscopy  at  a  depth  of  more 
than  three  millimeters(Shelton  and  Applegate  ,  Hu  et  al.  2010).  The  findings  imply  that  PAT 
has  potential  for  noninvasive  real  time  brain  mapping  of  cortical  processing  and  seizures. 

A  key  advantage  of  our  noninvasive  PAT  system  for  epileptiform  event  monitoring  is 
the  real  time  imaging  ability.  Seizures  are  generated  by  complex  interactions  of  a  large  group 
of  neurons  in  which  the  population  of  neurons  involved  varies  largely  from  moment  to 
moment(Schwartz  and  Bonhoeffer  2001).  By  achieving  1,000  times  faster  data  acquisition  of 
the  current  PAT  imaging  system  relative  to  our  previous  PAT  system  (Zhang  et  al.  2008),  we 
were  able  to  observe  that  the  spatial  aspects  of  the  seizure  focus  varied  over  time  as 

suggested  by  the  changes  in  hemodynamics  (Fig.  3c).  At  present,  the  imaging  speed  was 
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limited  only  by  the  10  Hz  laser  pulse-repetition-rate.  Since  considerably  faster  lasers  of  -500 
Hz  pulse-repetition-rates  are  now  commercially  available,  the  temporal  resolution  of  PAT 
will  soon  reach  a  level  of  a  few  milliseconds.  This,  coupled  with  the  compactness  of  an 
optical  system  will  make  it  possible  to  bring  a  PAT  system  to  the  bedside  for  real-time  and 
chronic  monitoring  of  seizures  in  infants  or  during  extra-operative  epilepsy  brain  mapping, 
for  example. 

One  hypothesis  is  that  long-standing  epilepsy  may  lead  to  the  development  of 
secondary  epileptogenic  regions  located  at  a  site  distant  from  the  original  focus,  a  factor  that 
may  reduce  the  likelihood  of  successful  epilepsy  presurgical  planning  (Cendes  et  al.  1995). 
The  ability  to  identify  the  epileptic  focus  and  associated  network  may  lead  to  better  epilepsy 
surgery  outcomes  for  many  individuals.  Various  methods  such  as  cross  correlation  and 
Granger  causality  have  been  explored  for  assessing  the  dynamic  directional  relationships 
among  brain  regions  using  time  series  data(Kaminski  et  al.  2001,  Bressler  et  al.  2008).  We 
used  Granger  causality  as  a  way  to  quantify  the  dynamic  interactions  amongst  the  primary 
seizure  focus,  secondary  seizure  foci,  and  middle  cerebral  artery  based  on  the  time-series 
photoacoustic  data.  Specifically,  Granger  causality  allowed  for  assessment  of  the  magnitude 
and  direction  of  temporal  relationships  during  overlapping  PAT  time-series  windows.  Results 
further  exemplified  the  temporal  aspects  of  seizure  onset,  seizure  secondary  spread,  and 
seizure  termination.  Collectively,  PAT  image  analysis  and  tools  for  effective  connectivity 
may  result  in  a  better  understanding  of  epileptic  networks  at  high  spatiotemporal  resolution. 

Little  data  are  available  regarding  a  link  between  vascular  changes  and  seizures 
(Ndode-Ekane  et  al.).  We  found  that  the  vasomotor  phenomena  of  micro  blood  vessels 
correlate  with  interictal  spike  and  wave  discharges  (Fig.  5d).  Interictal  spikes  generated  a 
strong  cerebral  metabolic  change  that  induced  cerebral  vessel  dilation  and  a  focal  increase  in 

cerebral  blood  flow  provided  oxygenated  hemoglobin  to  hypermetabolic  neurons,  with  a 
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corresponding  increase  in  total  hemoglobinfMa  et  al.  2009).  Vasomotion  during  seizure  onset 
and  spread  was  clearly  caused  by  the  oscillation  in  vessel  diameter  and  the  expansion  of  the 
vessel  cross  section.  We  reasoned  that  the  dilation  resulted  in  an  active  cerebral 
microvasculature  response  to  increased  metabolic  activity  accompanying  the  seizure(Myers 
and  Intaglietta  1976).  The  combination  of  multi-scale  imaging  ability  and  endogenous 
hemoglobin  contrast  makes  PAT  a  promising  tool  for  imaging  micro  vasculature  during 
seizure  onset.  Future  quantitative  PAT  studies  will  allow  us  to  resolve  many  of  the 
hemodynamic  components,  including  changes  in  level  of  hemoglobin  oxygenation  (i.e., 
deoxy-  and  oxy-hemoglobin),  cerebral  blood  flow,  and  the  rate  of  cerebral  oxygen 
metabolism,  for  a  detailed  understanding  of  the  neurovascular  coupling  phenomenon,  both 
during  and  in  between  seizures(Yuan  et  al.  2007,  Yao  et  al.  2009). 

The  present  work  represents  a  major  technological  and  scientific  breakthrough  in 
epilepsy.  Previous  attempts  to  track  from  moment  to  moment  populations  of  neurons 
participating  in  an  epileptiform  event  had  not  been  possible  with  any  technique  or  tool 
because  of  spatial  and/or  temporal  sampling  limitations.  The  high  spatial  and  temporal 
sampling  of  the  novel  PAT  system  allowed  for  the  first  time  the  complete  mapping  of  an 
epileptiform  event  in  vivo.  Another  major  breakthrough  is  that  this  is  the  first  report  of 
mapping  an  epileptiform  event  at  depths  well  below  the  cortex.  In  terms  of  impact,  this  is  the 
first  demonstration  of  emerging  optical  mapping  tool  in  the  surgical  evaluation  of  focal 
cortical  epilepsy,  since  accurate  localization  of  the  epileptic  focus,  propagation  paths,  and 
subcortical  networks  critically  depends  on  the  precise  localization  of  epileptogenic  neurons 
and  networks.  Clearly  the  use  of  our  PAT  system,  experimental  paradigm,  results,  and 
analyses  advances  our  understanding  of  epilepsy  and  seizure  temporal  spatial  properties 
relative  to  previous  attempts.  The  noninvasive  yet  whole  surface  and  depth  capabilities  of 

the  PAT  system  allowed  for  the  first  time  to  actually  see  what  is  happening  during 
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ictogenesis  in  terms  of  seizure  onset  and  spread.  The  challenge  of  mapping  focal  epilepsy 
stems  from  the  observation  that  the  pathology  associated  with  focal  epilepsy  is  often 
distributed  across  a  number  of  brain  sites.  Seizures  in  animal  models  and  in  people  often  have 
a  multifocal  or  broadly  synchronized  onset.  The  implication  of  these  observations  is  that  each 
of  the  sites  could  act  independently  to  initiate  a  seizure  or,  potentially,  to  drive  another  site 
into  a  seizure.  The  current  study  lends  support  to  the  theory  that  seizure  onset  and  spread 
involves  a  rich  interplay  between  multiple  cortical  and  subcortical  foci  during  the  onset  and 
spread  of  focal  epilepsy.  The  findings  are  timely  in  the  sense  that  the  neuroscience 
community  is  questioning  the  long-held  dogma  that  seizure  onset  involves  some  sort  of 
single  focus  or  epicenter  in  favor  of  the  emerging  thought  that  seizures  instead  involve 
multiple  cortical  and  subcortical  foci. 
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Figure  Legends 

Figure  1  |  Real-time  PAT  system  for  seizure  dynamics,  (a.)  Schematic  of  the  real  time 
PAT  system.  A  192-element  full-ring  transducer  array  was  used  to  capture  the  PA  signal 
during  seizure  onset.  Using  64-channel  parallel  data  acquisition  and  3:1  electronic 
multiplexing,  imaging  speed  was  up  to  0.33sec/frame  (Supplemental  Figure  1).  The  spatial 
resolution  of  this  imaging  system  was  — 150  microns  (Supplemental  Figure  2).  (b.) 
Noninvasive  PAT  imaging  of  a  rat  brain  in  vivo  with  the  skin  and  skull  intact.  MCA,  middle 
cerebral  artery;  RH,  right  hemispheres;  LH,  left  hemispheres;  LOB,  left  olfactory  bulbs; 
ROB,  Right  olfactory  bulbs,  (c.)  Open-skull  photograph  of  the  rat  brain  surface  acquired  after 
the  PAT  experiment.  Numbers  1-5  indicate  the  corresponding  blood  vessels  in  the  PAT 
image  and  rat  brain  photograph. 

Figure  2  |  Noninvasive  epileptic  foci  localization,  (a.)  EEG  recordings  of  seizure  onset  after 
BMI  injection,  (b.)  Schematic  showing  the  location  of  BMI  injection  and  saline  injection 
(control).  B,  BMI  injection;  S,  saline  injection,  (c.)  Reconstructed  PAT  image  before  (left) 
and  after  (right)  the  BMI  injection.  The  BMI  and  saline  were  injected  into  the  left  and  right 
parietal  neocortex,  respectively.  A  significant  increase  of  optical  absorption  is  seen  in  the 
region  of  the  BMI  injection,  while  no  absorption  contrast  is  observed  in  the  region  of  the 
saline  injection  relative  to  its  surroundings,  (d.)  Series  of  PAT  images  from  three 
representative  transverse  planes  parallel  to  the  skin  surface  with  755  nm  wavelength.  Slices 
5,  7,  and  9  are  the  images  obtained  3,  5  and  7  mm  under  the  skin,  respectively.  The  images 
show  different  spatial  patterns  at  the  seizure  foci  at  different  tomographic  layers,  (e.)  Three- 
dimensional  rendering  of  the  epileptic  foci  from  different  tomographic  layers  (Web 
Supplement,  Movie  SI). 
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Figure  3  |  Real-time  monitoring  of  ictal  onset,  (a.)  A  PAT  image  showing  the  cerebral 
vasculature  and  the  position  of  BMI  injection.  ROI,  region  of  interest;  Scale  bar:  2mm.  (b.) 
EEG  recordings  showing  the  seizure  onset  about  1  min  after  BMI  injection,  (c.)  PAT  images 
recording  the  ictal  onset  in  real  time.  At  lmin  6.267  s  the  area  of  the  focus  was  0.2053  mm2. 
The  area  of  the  focus  then  increased  in  size  over  the  next  2  s  to  0.5004mm2,  corresponding  to 
an  increase  in  the  amplitude  of  the  EEG  spikes.  The  area  of  the  seizure  focus  was  derived 
from  the  PAT  images  by  thresholding  to  a  pixel  value  one  standard  deviation  above  the  pixel 
values  from  the  area  of  the  focus  during  the  control  conditions.  Scale  bar:  500pm  (Web 
Supplement,  Movie  S2). 

Figure  4  |  In  vivo  Maps  of  mirror  foci  propagation,  (a.)  Left,  the  PA  image  of  the  rat  brain; 
the  yellow  dotted  rectangular  shows  the  ROI  for  analysis;  Right,  noninvasive  imaging  of  the 
primary  and  secondary  seizure  foci  in  bilateral  homotopic  cortex.  ai_4  shows  how  mirror  foci 
were  generated  and  diminished.  The  PA  signal  in  the  secondary  seizure  foci  was  smaller  in 
magnitude  and  delayed  in  time  compared  to  the  signal  recorded  from  the  primary  foci,  (b.) 
Left,  Granger  analysis  of  the  image  around  the  foci  and  its  ‘surround’  area;  Right,  PA  signal 
from  the  foci  (red)  and  surround  (blue)  during  an  ictal  event.  The  PA  signal  surrounding  the 
foci  (blue  line)  shows  inverted  signal  compared  to  that  from  the  seizure  foci  (red  line) 
corresponding  to  neuronal  inhibition,  (c.)  Granger  analysis  of  the  seizure  foci,  secondary 
seizure  foci  and  middle  cerebral  artery,  ci,  a  PAT  image  showing  the  nine  selected  ROIs;  Ci, 
connectivity  relation  between  the  primary  and  secondary  seizure  foci;  C3,  connectivity 
relation  between  the  primary  focus,  secondary  seizure  foci  and  middle  cerebral  artery.  The 
line  width  stands  for  the  relative  Granger  causal  influence  strength. 

Figure  5  |  Changes  of  blood  vessel  diameter  during  interictal  onset,  (a.)  PA  image  of  the 
rat  cortex  vasculature  with  the  intact  scalp  and  skull.  The  red  arrow  indicates  the 

microvascularity  (MV)  along  the  yellow  dashed  line  direction  for  a  cross-sectional  scanning 
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using  a  50  MHz  ultrasound  transducer,  (b.)  Typical  photoacoustic  signals  from  the  targeted 
vessel  at  different  times,  (c.)  EEG  recordings  show  the  interictal  spikes  and  (d.)  Change  of 
the  vessel  size  captured  by  PAT.  There  was  a  clear  correlation  between  the  interictal  spikes 
and  the  changes  of  the  vessel  size. 

Supplement  Figures 

Figure  1  |  Temporal  resolution  demonstration  of  the  PAT  system.  Photoacoustic  images 
depicting  frames  from  a  10-second  sequence  of  dynamic  ink  flow  through  a  0.3  mm  diameter 
tube  (Web  Supplement,  Movie  S3).  The  movie  shows  that  the  temporal  resolution  of  this 
imaging  system  was  -0.33  seconds/  frame. 

Figure  2  |  Spatial  resolution  measurement  of  the  PAT  system,  (a.)  Cross-sectional  PA 
image  of  a  two-copper-rods  phantom.  The  two  copper  rods  (diameter:  0.05mm)  were 
embedded  in  a  cylindrical  tissue-mimicking  phantom  at  a  depth  of  1 0  mm.  The  center-to- 
center  distance  between  the  two  copper  rods  was  approximately  2.2  mm.  (b.)  The  normalized 
profile  of  the  reconstructed  PA  image  shown  in  (a.)  along  y=10  mm.  The  half-  and  quarter- 
amplitude  widths  were  obtained  by  measuring  the  distance  between  points  AB  and  between 
CD,  respectively.  The  spatial  resolution  of  the  imaging  system  calculated  with  Rayleigh’s 
law  was  150pm. 
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